# A multi-state model incorporating estimation of excess hazards and multiple
time scales

Caroline E. Weibull  Department of Medical Epidemiology and Biostatistics,
Karolinska Institutet, Box 281, 171 77 Stockholm, Sweden  Department of
Medicine Solna, Clinical Epidemiology Division, Karolinska Institutet,
Eugeniahemmet, T2, Karolinska Universitetssjukhuset, 171 76 Solna, Sweden
Paul C. Lambert  Department of Medical Epidemiology and Biostatistics,
Karolinska Institutet, Box 281, 171 77 Stockholm, Sweden  Biostatistics
Research Group, Department of Health Sciences, University of Leicester,
University Road, Leicester, LE1 7RH, UK  Sandra Eloranta  Department of
Medicine Solna, Clinical Epidemiology Division, Karolinska Institutet,
Eugeniahemmet, T2, Karolinska Universitetssjukhuset, 171 76 Solna, Sweden
Therese M.L. Andersson  Department of Medical Epidemiology and Biostatistics,
Karolinska Institutet, Box 281, 171 77 Stockholm, Sweden  Paul W. Dickman
Department of Medical Epidemiology and Biostatistics, Karolinska Institutet,
Box 281, 171 77 Stockholm, Sweden  Michael J. Crowther  Department of Medical
Epidemiology and Biostatistics, Karolinska Institutet, Box 281, 171 77
Stockholm, Sweden  Biostatistics Research Group, Department of Health
Sciences, University of Leicester, University Road, Leicester, LE1 7RH, UK

###### Abstract

As cancer patient survival improves, late effects from treatment are becoming
the next clinical challenge. Chemotherapy and radiotherapy, for example,
potentially increase the risk of both morbidity and mortality from second
malignancies and cardiovascular disease. To provide clinically relevant
population-level measures of late effects, it is of importance to (1)
simultaneously estimate the risks of both morbidity and mortality, (2)
partition these risks into the component expected in the absence of cancer and
the component due to the cancer and its treatment, and (3) incorporate the
multiple time scales of attained age, calendar time, and time since diagnosis.
Multi-state models provide a framework for simultaneously studying morbidity
and mortality, but do not solve the problem of partitioning the risks.
However, this partitioning can be achieved by applying a relative survival
framework, by allowing is to directly quantify the excess risk. This paper
proposes a combination of these two frameworks, providing one approach to
address (1)-(3). Using recently developed methods in multi-state modeling, we
incorporate estimation of excess hazards into a multi-state model. Both
intermediate and absorbing state risks can be partitioned and different
transitions are allowed to have different and/or multiple time scales. We
illustrate our approach using data on Hodgkin lymphoma patients and excess
risk of diseases of the circulatory system, and provide user-friendly Stata
software with accompanying example code.

##  1 Introduction

Marked improvement in cancer patient survival in recent decades has increased
the importance of studying treatment side-effects[1], i.e., late effects. For
example, chemotherapy and radiotherapy increase the risk of both secondary
malignancies and cardiovascular disease (CVD)[2, 3]. As therapies are
developed and adapted to maximise the chances of cure while reducing risk of
treatment toxicity, it is of interest to simultaneously estimate real-world
morbidity and mortality. Today, a major challenge in the field of late effects
lies in estimating the proportion of the observed morbidity and mortality that
can be ascribed to the underlying cancer and its treatment. Hence, the
standard approach is to simply estimate the cumulative observed incidence and
mortality[4]. This article proposes a population-based approach to this
problem bringing together statistical methods for relative survival and multi-
state modelling. The former allows estimation of excess incidence, which under
certain assumptions corresponds to the incidence above and beyond that
expected in the absence of cancer and treatment thereof, and the latter
enables investigation of patient trajectories going from an initial state
(diagnosis), potentially via intermediate states (e.g., non-fatal illnesses
such as CVD), to final absorbing states (different causes of death). Our focus
is on the illness-death model where the illness state can be experienced also
in the absence of cancer (unlike recurrence, which only occur in cancer
patients). In addition, we know only when patients enter each state and not
whether it was due to the cancer. We are interested in partitioning the
probability of being in the illness state into that due to the cancer and its
treatment, and that which would be expected in the absence of cancer.

Previous attempts to combine relative survival with multi-state modelling has
proposed a Markov piecewise constant rate model incorporating estimation of
excess mortality[5, 6]. This method has been further generalised by Gillaizeau
et al. to include semi-Markov models[7]. In these papers, the relative
survival methods are applied to absorbing states alone, i.e., excess mortality
can be estimated but not excess incidence. Moreover, focus has been on
relative effects rather than predictions. Our suggested approach enables
partitioning of both excess morbidity and mortality in a Markov or semi-Markov
setting, and we are unaware of existing methods for doing this.

In population-based cancer patient survival, relative survival is a commonly
used summary measure of patient survival as it captures both the direct and
indirect contribution of a cancer diagnosis on mortality by comparing the
observed survival of the patients to the expected survival in a comparable
cancer-free population[8, 9, 10]. In the standard setting, the outcome is all-
cause death and the mortality rates expected in the absence of cancer are
typically retrieved from publicly available population mortality tables.
However, the same approach can be used to estimate excess incidence rates[11].
As such, the excess incidence of some disease of interest can be interpreted
as the portion of incidence attributable to the underlying cancer and its
treatment, given that the expected rates are appropriate for the population
diagnosed with cancer.

Lambert et al. have previously developed methods for estimating cumulative
probabilities of death parametrically in a relative survival framework [12].
Eloranta et al. extended these methods to enable partitioning of the excess
mortality associated with cancer into component parts such as excess
cardiovascular mortality and remaining excess cancer mortality[13]. This has
been applied to study temporal trends in excess mortality due to diseases of
the circulatory system (DCS) among Hodgkin lymphoma (HL) patients[14]. Further
extensions of these methods have been applied to study excess DCS incidence in
HL patients[11]. In the latter paper, the expected incidence rate was modelled
using a constant rate model with two time scales (calendar year and age) with
twoway interactions. In the mentioned papers, the cumulative probabilities
were calculated using numerical integration, in the presence of competing
risks. However, considering morbidity as an absorbing state, as in the
application on HL and excess DCS incidence, prevents identification of fatal
complications, which in turn makes comparisons with previous studies on excess
mortality more difficult.

Multi-state models are a generalisation of competing risks models and one way
to study patient trajectories[15, 16]. The most commonly used method to
estimate transition rates in a multi-state model is using a semi-parametric
approach, such as the Cox model. In a recently published study, Crowther and
Lambert[17] suggested a framework for transition-specific parametric
distributions in multi-state models. In general, parametric survival models
have several advantages over semi-parametric models, including the potential
of modelling more than one timescale, more efficient incorporation of time-
dependent effects, ease of prediction (both in- and out-of-sample), and
greater flexibility in quantifying absolute risks. In addition, estimation of
uncertainty is easier, and they enable extrapolation of survival over the full
life span (although this should always be done with caution). In the
previously proposed method[17], a separate parametric model can be fitted to
each transition, which allows the use of transition-specific distributions.
The transition probabilities are calculated using a simulation approach[18,
19]. Due to the flexibility of the simulation approach, other quantities of
interest can be easily predicted, such as the average time spent in each
state.

In survival analysis, there may be more than one plausible timescale on which
to quantify an event rate[20]. An obvious example is how incidence,
prevalence, and mortality due to various cancers vary as functions of attained
age and calendar time. Acknowledging this may introduce transition-specific
and/or multiple timescales, which are rarely investigated since the
computational complexity increases substantially. Work in this area has been
sparse, particularly moving to competing risks and the more general multi-
state setting. Lee and Fine[21] developed cause-specific Cox models in a
competing risks setting, allowing cause-specific timescales, where attained
age was the underlying timescale for one cause of death and time since
diagnosis was assumed for the competing cause of death. Wolkewitz et al.[22]
developed a two timescale Cox model applied to modelling risk of infection
among patients admitted to Intensive Care Units (ICU), with discharge and
death as competing risks, using time since entry to the ICU and calendar time
as timescales. In a parametric framework, the most substantial contribution to
the field is by Iacobelli and Carstensen[16], who provided a general approach
for estimating multi-state models with multiple timescales, centred on time-
splitting and estimating rates within a parametric spline-based modelling
framework using Poisson regression. Predictions of transition probabilities
were obtained using simulation, and confidence intervals calculated using a
parametric bootstrap.

Motivated by the need to simultaneously investigate excess morbidity and
mortality among cancer survivors, this study proposes one solution to several
methodological challenges related to the clinical application. Starting from
an illness-death model (with cancer diagnosis as the initial starting state)
we incorporate estimation of excess hazard rates into the transition model for
the illness state. Moreover, in the expected hazard rate model we allow for
multiple time scales (calendar time and age) while remaining transitions have
one single time scale (time since study entry or time since entry to current
state). Generalisation beyond the illness-death model is easily achieved. The
approach is illustrated on a study of excess incidence and mortality from
diseases of the circulatory system (DCS) among patients treated for HL.

The article is structured as follows. In Section 2 we give a brief overview of
recently developed methods for multi-state modelling in the cause-specific
framework, after which we describe our proposed approach for each transition.
Section 3 provides an illustrative example on patients diagnosed with HL and
excess DCS morbidity and mortality. We conclude with a discussion and suggest
possible extensions in Section 4.

##  2 Multi-state modelling incorporating estimation of excess hazards

We begin by conceptually describing the method proposed by Crowther and
Lambert[17]. We define a global vector of available covariates, $\bm{X}$,
since the same covariates may be used in multiple models. The fitted
transition-specific models can be as complex as required (with respect to
functional form of covariates, interactions, and time-varying effects) and one
can use different parametric models for different transitions (e.g., a
flexible parametric for one transition and a Weibull model for another). Once
the set of transition models has been chosen, transition probabilities are
calculated for a covariate pattern of interest. The details have been
extensively described elsewhere[17] but in short, point estimates are obtained
by simulating a large sample of patient trajectories ($n$) through the fitted
multi-state model. Confidence intervals are calculated from repeated
simulations ($m$), each drawing from a multivariate normal distribution with
the estimated parameter vector as mean and the variance-covariance matrix (of
the parameters) as variance, i.e., a parametric bootstrap.

Time scales: $t$: Time since cancer $a$: Age at cancer diagnosis $c$: Year of
cancer diagnosis $u$: Time of illness State 1: Alive with cancer State 2:
Expected illness Excess illness State 3: Dead before illness State 4: Dead
after illness $h_{1}^{*}(t+a,t+c)$ $\lambda_{2}(t)$ $h_{3}(t)$ $h_{4}(t-u)$
Figure 1: Illustration of suggested approach. $h_{1}^{*}(t+a,t+c)$ refers to
the illness incidence rate that is expected in the absence of cancer with
attained age ($t+a$) and calendar year ($t+c$) as time scales,
$\lambda_{2}(t)$ is the excess incidence of illness with time since cancer
diagnosis ($t$) as the time scale, $h_{3}(t)$ is the mortality rate with time
since diagnosis ($t$) as the time scale, and $h_{4}(t-u)$ is the mortality
rate among cancer survivors who experience the illness, with time since
illness ($t-u$) as the time scale.

Figure 1 shows a representation of our suggested extension. While a standard
illness-death model consists of three states (healthy, ill, dead), we
partition the illness state into two unobserved states - excess and expected
illness (dashed lines). We further have two death states, death following
illness, and death without having experienced the illness, which allows us to
quantify the impact of the illness under study on mortality. In total, we have
four states, of which one is partitioned into two components. Patients start
out in state 1; alive after a diagnosis of cancer. Patients might then
experience some illness (state 2), which can be either related to the cancer
and its treatment (“excess”) or not (“expected”). The transition rate to the
expected illness component of state 2, $h_{1}^{*}(t+a,t+c)$, is modelled using
population records of disease incidence with attained age ($t+a$) and calendar
time ($t+c$) as time scales. The transition to the excess illness component of
state 2, $\lambda_{2}(t)$, is estimated using a parametric excess hazards
model, incorporating the modelled expected illness rate. Patients might die
prior to experiencing the illness under study (state 3). This mortality rate,
$h_{3}(t)$, can be modelled as all-cause or partitioned into cause-specific
mortality models. The time scale for both the excess incidence model and this
mortality model is time since cancer ($t$). Some of the patients experiencing
the illness (excess or expected) will transition to the absorbing state of
death (state 4). This transition, $h_{4}(t-u)$, can be modelled as an all-
cause or cause-specific mortality model, with either one time scale (here we
illustrate using time since illness, $t-u$) or multiple time scales (a natural
extension would be to add $t$). Any of the models used for transitions into
states 3 and 4 can be extended to excess hazards models, given the appropriate
expected rate information.

###  2.1 Transition 1: Modeling the expected incidence

In standard relative survival we use mortality rates from publicly available
life tables (usually stratified by year, age, and sex) as a measure of the
expected mortality in the general population. In order to incorporate such an
expected rate into the transition-specific model framework, and hence the
simulation approach used to obtain predictions, we must model the rate with an
appropriate parametric specification. As incidence, rather than mortality, is
the entity of interest here, population incidence tables are required to model
the expected incidence rate, i.e., the number of events and person-time for
each combination of the stratification variables, as shown in Table 1, are
needed.

    yearΨ  sex      age       d         y
    2001     1       18      62   74115.36
    2001     1       19      83   74258.12
    2001     1       20      85   75230.58
    ...
    2002     1       18      63   70099.6
    2002     1       19      64   73899.11
    ...

Table 1: Example of data setup for modelling the expected incidence rate of
illness. For every combination of calendar year, sex, and age, the number of
events (d) and person-time (y) are calculated.

While any parametric hazards model can be used to estimate the (log) expected
rate, we illustrate how to model the transition from state 1 to the expected
illness component of state 2 using a spline-based Poisson regression model.

$\log\left(\frac{d}{y}\right)=f(a)+g(c)+\mathbf{x}_{1}\mathbf{\beta}_{1}$ (1)

As disease incidence in a population is typically a function of both age and
calendar year, both are used as time scales here. Additionally, population
incidence/mortality tables are typically not available on individual level,
which precludes modelling the rate with one of these variables as the time
scale and other as a time-fixed covariate. Age and calendar year can be
modelled using a range of functions $f(\cdot)$ and $g(\cdot)$, for example,
assuming linearity or using restricted cubic splines or B-splines. In this
paper, we use restricted cubic splines to model both timescales, which are
easily conducted within a Poisson framework[16].

We will use $\mathbf{x}_{1}$, where $\mathbf{x}_{1}\subset\mathbf{X}$, to
denote the vector that contains additional stratification factors for the
population incidence table; this typically only contains sex but might also
include race, region, or social class. We further let $\mathbf{\beta}_{1}$
denote the vector of associated coefficients. For simplicity, we assume
proportional hazards, although this assumption can be relaxed by including
interaction terms between the time scales and covariates into the model.

Equation 1 specifies the model for the log expected rate as a function of the
two continuous time scales and any additional covariates. However, for
implementation purposes, the model needs to be specified in terms of the main
time scale for the multi-state model (time since cancer diagnosis):

$\log
h^{*}_{1}(t+a,t+c|\mathbf{x}_{1})=f(t+a)+g(t+c)+\mathbf{x}_{1}\mathbf{\beta}_{1}$
(2)

We provide more details on the implementation in the Appendix.

###  2.2 Transition 2: Modeling the excess incidence

While the expected incidence rate is modelled using population incidence data
on a grouped level, the remaining models use the patient cohort data. The
setup can be either in wide format (each patient has one row with information
on all transitions) or long format as described by Putter et al.[15] (each
patient has a separate row of data for every transition). Long format (i.e.,
stacked data notation) can be useful in situations where covariate effects are
shared across several transitions. We model the excess incidence of illness
using a parametric relative survival model. The most appropriate choice of
time scale for the excess incidence is time since cancer diagnosis ($t$).
Using the notation from Figure 1, the overall incidence, $h(t,t+a,t+c)$, among
the patients can be written as the sum of the incidence expected in the
absence of cancer, $h^{*}_{1}(t+a,t+c)$, and the excess incidence,
$\lambda_{2}(t)$, associated with the cancer and its treatment:

$h(t,t+a,t+c|\mathbf{x}_{1},\mathbf{x}_{2})=h^{*}_{1}(t+a,t+c|\mathbf{x}_{1})+\lambda_{2}(t|\mathbf{x}_{2})$
(3)

The vector $\mathbf{x}_{2}$, where $\mathbf{x}_{2}\subset\mathbf{X}$,
represents the set of covariates related to the excess rate, which usually
includes sex, baseline age and calendar year. We choose to model the
transition from state 1 to the excess illness component of state 2 on the log
cumulative excess hazard scale. The overall cumulative hazard,
$H(t,t+a,t+c|\mathbf{x}_{1},\mathbf{x}_{2})$, can be expressed as the sum of
the cumulative expected hazard, $H^{*}_{1}(t+a,t+c|\mathbf{x}_{1})$, and the
cumulative excess hazard, $\Lambda_{2}(t|\mathbf{x}_{2})$:

$H(t,t+a,t+c|\mathbf{x}_{1},\mathbf{x}_{2})=H^{*}_{1}(t+a,t+c|\mathbf{x}_{1})+\Lambda_{2}(t|\mathbf{x}_{2})$
(4)

Although any parametric relative survival model can be used, we will
illustrate using a flexible parametric relative survival model[23], for which
the cumulative excess hazard is defined as:

$\log\Lambda_{2}(t|\mathbf{x}_{2})=\log\Lambda_{2,0}(t)+\mathbf{x}_{2}\mathbf{\beta}_{2}.$
(5)

For illustrative purposes, we choose to model with only time since cancer
diagnosis as the time scale and we incorporate age and year of diagnosis as
baseline covariates (included in $\mathbf{x}_{2}$), but it is possible to
incorporate these as secondary time scales. Moreover, interactions, including
time-dependent effects, can easily be included. Given the completeness of
incident cases of the illness under consideration, the predicted expected rate
is treated as fixed, i.e., assuming no uncertainty, tailored at the individual
patient level and incorporated into the likelihoood [24].

###  2.3 Transition 3: Modeling the mortality rate (before illness)

In addition to being at risk for late complications, patients are at risk of
dying from the cancer itself and from other causes unrelated to the cancer. In
the simplest case, interest is not in the mortality from these causes (other
than to account for this as a competing event) and any standard parametric
all-cause mortality model can be used to model the transition from state 1 to
state 3. Note however, that this will implicitly not include deaths caused by
the illness under study. We illustrate using a flexible parametric survival
model[25], where modelling is done on the cumulative hazard scale:

$\log H_{3}(t|\mathbf{x}_{3})=\log
H_{3,0}(t)+\mathbf{x}_{3}\mathbf{\beta}_{3}$ (6)

The log cumulative baseline hazard, $\log H_{3,0}(t)$, is modelled using
splines, much like in the excess incidence model above, and we have our
covariate vector $\mathbf{x}_{3}\subset\mathbf{X}$ with associated
coefficients, $\mathbf{\beta}_{3}$. Although mortality due to other causes
would be best modelled using attained age as the time scale, the mortality
among these patients is highly driven by cancer deaths, and thus time since
cancer diagnosis is the more reasonable choice of time scale. However, we
might also extend this model to incorporate multiple time scales. Though not
described in detail here, extending this to a more complex model - such as an
excess mortality model - is possible. This would introduce also a model for
the expected mortality from all causes of death except for the illness.

###  2.4 Transition 4: Modeling the mortality rate (after illness)

Once patients enter the illness state, it is not possible to separate between
patients who experienced the illness due to the cancer and patients who would
have been diagnosed with the illness regardless of the cancer. So even though
there are two transitions going into state 2, i.e., into the expected and
excess illness components of the illness state, there is just one transition
going from state 2 to whichever states patients can experience next. As with
the other transitions, any parametric model can be applied. The choice of time
scale can be either time since cancer diagnosis (resulting in a delayed entry
model) or time since illness diagnosis, i.e. a clock reset model. It is also
possible to use both as time scales. We again illustrate with a flexible
parametric survival model, now with time since illness as the time scale.

$\log H_{4}(t-u|\mathbf{x}_{4})=\log
H_{4,0}(t-u)+\mathbf{x}_{4}\mathbf{\beta}_{4}$ (7)

The log cumulative baseline hazard, $\log H_{4,0}(t-u)$, is modelled using
splines, and we have our covariate vector $\mathbf{x}_{4}\subset\mathbf{X}$
with associated coefficients, $\mathbf{\beta}_{4}$. Similar to the transition
to death before illness, a more complex model can be incorporated. When
applying an excess rate model, the excess mortality could be further
partitioned into excess illness- and cancer mortality following the method
suggested by Eloranta et al. (2012)[13]. However, this can become complex as
it requires calculating expected mortality rates conditioning on the illness
event.

For implementation purposes, the excess and expected components that make up
state 2 must be treated as separate states (say 2a and 2b) when defining the
transition matrix, since they have separate transition hazards governing their
rate of occurrence. Subsequently, the transition rates for the next
transitions of going from 2a to 3, and 2b to 3, must share the same transition
model, as the partitioning of state 2 is essentially theoretical. As a result
of this, the accompanying software package has been extended to allow states
to share transition model.

##  3 Illustrative example

We illustrate the suggested approach with an application to the modelling of
excess incidence and mortality from diseases of the circulatory system (DCS)
among patients treated for HL. HL is often pointed out as a model disease in
terms of improvements in survival, going from a fatal disease in the mid-1900
to around 90% of young HL patients today being cured from their lymphoma[26,
27]. These improvements are primarily due to major advances in the management
of these patients, including more accurate radiotherapy, effective multi-agent
chemotherapy, and improved staging procedures. With more patients surviving
their HL, research focus has somewhat shifted towards understanding and
reducing treatment toxicity. The two major complications with potential fatal
outcome are secondary malignancies and DCS.

From the Swedish Cancer Register, we identified 4,479 patients diagnosed with
HL between 1985 and 2013, aged 18-80 at diagnosis. This cohort has been
described in detail previously[11]. Patients were followed until death,
emigration, or end of study period ($31^{\text{st}}$ of December, 2014),
whichever came first. Incident cases of DCS were linked to the HL patients
using the Swedish Inpatient Register. Only events with DCS as the main
diagnosis were considered. At the end of each individual’s follow-up, 2,556
patients had not experienced any of the secondary events, 978 patients died
due to any cause prior to experiencing DCS, and 945 patients had been
diagnosed with DCS. Among these, 425 were alive and living with DCS at the end
of study, and 520 patients had died (from any cause of death) (Figure 2).

HL Start=4,479 End=2,556 DCS Start=0 End=425 Dead before DCS Start=0 End=978
Dead after DCS Start=0 End=520 $n=945$$n=978$$n=520$ Figure 2: Illustration of
patient trajectories for the 4,479 patients diagnosed with HL in Sweden
between 1985 and 2013, ages 18-80 years at diagnosis. ‘Start’ represents the
number of patients who begin in each state, and ‘End’ represents the number of
patients remaining in each state at the end of follow-up.

For DCS, no population incidence table with sufficiently detailed data is
publicly available in Sweden. Instead, a cohort of approximately 10 million
individuals (Swedish residents who were either born after 1931 and residing in
Sweden after 1961, and their parents, or who were included in a census
1960–1990) was constructed using several population-based registers, as a
representation of the Swedish population. Data were stratified by calendar
year, age, and sex, and the number of DCS events and amount of person-years at
risk were calculated. The expected DCS rate was modelled using Equation 2,
with 5 degrees of freedom for each time scale and incorporating sex.

For modelling all remaining transitions, flexible parametric survival models
were used with year of diagnosis, age at diagnosis, and sex included as
baseline covariates. Year and age were parameterised using restricted cubic
splines with 5 degrees of freedom in the models for excess DCS incidence and
mortality before DCS. For post-DCS mortality, 3 degrees of freedom were used
due sparseness of data. No parameter effects were shared between models. For
simplicity, only main effects of covariates were considered, however,
interaction terms could easily be incorporated.

The transition from HL to excess DCS was modelled with time since HL diagnosis
as the time scale. The predicted expected rate was incorporated into the
likelihood to enable estimation of the excess rate. For the transition from HL
to dead before DCS, time since HL diagnosis was used as time scale and all-
cause mortality was the outcome of interest. Patients with DCS identified from
their cause of death certificate ($n=100$) were given that date as their
diagnosis date and the following day as their death date. Because of this, the
all-cause mortality before DCS implicitly did not include any deaths caused by
DCS. Death following a DCS diagnosis was modelled with time since DCS as the
underlying time scale (clock-reset).

Figure 3: Estimated transition rates for male patients aged 30 and 60 years at
diagnosis. Expected DCS incidence was modelled as a function of calendar time.
Remaining panels are for patients diagnosed in 1995, with excess DCS incidence
and other-cause (excluding DCS) mortality before DCS as functions of time
since cancer diagnosis, and all-cause mortality after DCS as a function of
time since DCS diagnosis.

Figure 3 shows the transition rates for male patients aged 30 and 60 year old
at diagnosis. The expected DCS rate is presented as a function of calendar
time, while the remaining rates are for patients diagnosed in 1995, and
plotted over time since diagnosis (cancer or DCS). The expected rates declined
slightly following the mid-1980s, and the rate among younger patients was very
low throughout. The excess DCS rate and both mortality rates were all high
early during follow-up.

Using the transition-specific model estimates, point estimate of the
transition probabilities were calculated by simulating 1,000,000 subjects
through the multi-state model. Figure 4 shows the stacked transition
probabilities over time since HL diagnosis for male patients diagnosed in 1995
at ages 30 (left panel) and 60 (right panel) years at diagnosis. For 30-year-
olds, the probability of remaining in the HL state at 15 years after cancer
diagnosis was approximately 80% whereas for 60-year-olds this probability was
just above 20%. The probability of being diagnosed with DCS within the first
15 years as expected in the absence of HL and its treatment, and still being
alive, was very small for younger patients. The probability of death (before
or after a DCS diagnosis) was expectedly much larger for older patients.

Figure 4: Probability of being in each state for male patients diagnosed in
1995, aged 30 and 60 years, as a function of time since HL diagnosis.

Although stacked graphs provide a useful tool for illustrating in which state
the patients are at different points in time, they lack information on
uncertainty. As no parameters were shared between transition models, re-
sampling to obtain confidence intervals was done using the transition-specific
coefficient vectors and associated variance-covariance matrices, with
$n=10,000$ and $m=1,000$. Figure 5 shows the transition probabilities with
associated confidence intervals for HL patients diagnosed in 1995 at ages 30
and 60 years.

(a) Age 30 years at HL diagnosis

(b) Age 60 years at HL diagnosis

Figure 5: Probability (with 95% confidence interval) of being in each state
for male patients diagnosed in 1995, aged 30 (a) and 60 (b) years, as a
function of time since diagnosis.

The fitted models can further be used to predict other clinically relevant
measures, such as differences in transition probabilities across varying
covariate patterns. Figure 6 illustrates the differences in transition
probabilities between females and males diagnosed in 1995, aged 30 (top panel)
and 60 (bottom panel) years.

(a) Age 30 years at HL diagnosis

(b) Age 60 years at HL diagnosis

Figure 6: Difference in probabilities (with 95% confidence interval) of being
in each state between female and male patients diagnosed in 1995, aged 30 (a)
and 60 (b) years, as a function of time since diagnosis. Abbreviations; M,
males; F, females.

Another measure of clinical importance is the proportion of excess DCS
(calculated as the probability of excess DCS over the probability of any DCS),
which allows us to identify trends in the contribution of excess DCS on total
DCS risk as a whole. Figure 7 shows the proportion of excess DCS at ten years
after HL, for male patients diagnosed in 1985 and 2000, as a function of age
at diagnosis. We could alternatively quantify such a proportion over the other
timescales, i.e. as a function of time since diagnosis, or indeed year at
diagnosis.

Figure 7: Proportion excess DCS among total DCS probability at 10 years after
diagnosis, among male patients diagnosed in 1985 and 2000, by age at
diagnosis.

Example code to conduct the previous analysis is included in the Appendix.
Extended predictions such as length of stay in each state, the probability of
ever visiting a state, standardised predictions and further contrasts are also
provided.

##  4 Discussion

Research on late effects of treatment among cancer survivors is becoming
increasingly important as more patients are cured from their disease. Studying
life span trajectories requires complex methods, that consider follow-up not
just until the first event of interest, but through a series of intermediate
states until reaching an absorbing state. Advanced multi-state methods for the
cause-specific setting have been developed and are frequently used in applied
research, but are unable to separate between the probability of entering an
intermediate state associated with a diagnosis of cancer, and the probability
of entering the same state not related to cancer. Utilising recently developed
methods in multi-state modelling, where a simulation approach is taken
together with parametric transition models, we have proposed one way to
incorporate estimation of excess morbidity and mortality. This allows for a
separation between the predicted excess illness probability and expected
illness probability. We have further introduced the possibility for different
transitions to have multiple time scales within a relative survival setting.

The developments provided in this paper were motivated by the study of excess
morbidity and mortality of DCS among HL patients. The excess DCS mortality has
been shown to decline in Sweden since the mid-1970s[14], and is no longer a
major cause of death among this group of patients. Still, the risk of being
diagnosed with excess DCS persists in the new millennium among HL patients in
Sweden[11]. Although a reduction in excess DCS mortality does not necessarily
imply a reduction in excess incidence, it remains difficult to draw
conclusions on the relationship between HL and its treatment and risk of DCS
together with fatality of DCS. By incorporating morbidity and mortality in the
same setting, and separating between the incidence attributable to the cancer
and that expected in the absence of cancer, we have illustrated one way to
overcome this issue. Our approach has been implemented in the merlin and
`multistate` packages in Stata, and code illustrating how to produce the
results presented in this paper, can be found in the Appendix.

For illustrative purposes, we chose to simplify the models in our example.
Only main effects models were considered, and no information on clinical
covariates were included. Therefore, caution should be taken in terms of
interpreting the results presented herein. It is important to point out that
more complex, and realistic, models could be fitted to each of the
transitions. This includes, but is not limited to, interactions between
covariates and relaxing the proportional hazards assumptions [28]. As with any
standard parametric modelling process, diagnostics and goodness-of-fit
inference based on Wald tests, likelihood ratio tests, information criteria
(AIC, BIC) can be applied and used to guide model selection for each of the
transitions. Lee et al. (2017)[21] recently highlighted that standard
Martingales must be adapted for the cause-specific timescale setting, and as
such, further research is required to adapt them to our more general multiple
timescale setting. Indeed, the example model specified for the excess
incidence rate (Equation 3) here only depended on baseline age and calendar
year, whereas multiple timescales could be incorporated also for this model.
Moreover, the final predicted probabilities rely on four different models so
in all applications it is important to assess the goodness of fit of each
model. We exemplify this in the Appendix by varying the degrees of freedom;
3-5 degrees of freedom for transition 2 and 3, and 3-4 degrees of freedom for
transition 4 (Supplementary Figures 8 and 9). In our illustrative example we
partitioned the illness state into an expected and excess component and
focused on one of several possible late effects. However, partitioning the
absorbing states can be done in a similar manner given the appropriate
expected mortality table, and additional illness states, e.g., secondary
malignancies can be added to the model.

In terms of quantifying model results, we have given examples of how to
predict transition probabilities, difference in probabilities, and proportion
of excess illness over total probability of illness. Our approach can be used
to predict other quantities, such as length of stay in each state (including
differences and ratios thereof between patients with different covariate
patterns). All of the described extensions can be implemented with the
associated software package.

Several assumptions are made when using a relative survival framework.
Patients are assumed to be exchangeable with the general population,
conditioned on the stratification variables. If this assumption is not
fulfilled, population incidence or mortality tables require additional
stratification. If this is not possible, information on a control population
can be used to adjust the expected rate [29]. Additionally, if we wish to
interpret the excess rate as treatment-related, then we must assume that the
cancer is not directly associated with illness.

###  4.1 Alternative methods

As in most research applications, there exist multiple routes to achieve the
same goal – this is also the case here. An alternative to our chosen approach
would have been to design a matched cohort study (i.e., cancer patients would
be considered exposed and matched to healthy (unexposed) comparators, on age,
sex, and year). This is feasible in countries where appropriate control
individual population data is available, however, this is more often not the
case, and hence reference population incidence/mortality tables must be
utilised. As our focus is on quantifying excess (additive) risk, if a matched
cohort approach is taken, additive models would need to be fitted, with the
natural model choice being the Aalen additive model [30]. Although a
simplistic version of this may have been possible, that modelling approach has
substantial limitations which we believe strongly precludes its use. We are
strong advocates of flexible, yet parametric, modelling approaches. There is
limited research into the development of, for example, a spline-based Aalen
model, and even more so when attempting to incorporate (arguably essential)
time-dependent effects. Taking such an approach would require modelling a
baseline rate with additive multiple timescales (which, we would argue, should
be modelled in a multiplicative context), combined with an additional additive
linear predictor, would raise further challenges. As such, this led to the
proposed approach, where we model the baseline (expected) rate appropriately,
and combine it with an additive excess component by adopting a relative
survival approach. Adopting a flexible parametric approach is important to
then obtain smooth estimates of interpretable measures of risk from a
comprehensive multi-state model. Having said that, developing an analysis
approach for a matched cohort design with multiplicative multiple timescales,
and an additive excess component is an important topic for further research,
which we are now actively developing as a follow-on study to this paper.

###  4.2 Limitations and future research

In our example, we only considered patients’ first occurrence of DCS. In
reality, patients may experience multiple DCS occurrences prior to death, or
indeed other intermediate events of interest. Incorporating recurrent DCS
events raises further challenges, for example, a first occurrence of a DCS
event could be cancer-related, whereas a second is not. By introducing
subsequent intermediate events into the model framework (with potential
further partitioning into excess and expected components), we can accommodate
such a setting. This extension is currently being considered.

One could argue that if population level predictions are the entity of
interest, marginal predictions rather than covariate-specific, as we have
presented, should be derived. Such marginal predictions involve computing the
statistic of interest at all observations’ covariate patterns, and taking the
average. We have implemented them within the associated software through the
standardize option, which is highly computationally intensive in the current
implementation, and as such requires further development.

We employed a latent time approach for our simulation method of prediction,
also used in Crowther and Lambert (2017)[17]; however, Bluhmki er al.
(2019)[31] have highlighted the disconnect between simulating multiple
unobservable event times and the real world setting of only observing one
outcome. We agree, but would argue that such a hypothetical construct does not
preclude the use of the algorithm, and as pointed out by Allignol et al.
(2011)[32], quantities obtained from either method are computationally valid.
Both methods are now available within the associated software. To derive
confidence intervals for all quantities presented in this paper, we have
utilized the parametric bootstrap. Alternative techniques, such as the non-
parametric bootstrap, may be beneficial as it is robust to model
misspecification [33]. However, this would be more computationally intensive
due to the necessity of re-fitting all transition models on each bootstrap
sample. Finally, we may derive standard errors utilising the multivariate
delta method to further avoid computationally intensive bootstrapping. Such
investigations deserve further research.

##  5 Conclusion

Cancer survivors are at risk of a range of diseases for the remainder of their
lives. To support this clinical challenge with scientific real-world evidence,
we believe that our approach can be useful for predicting survivors
trajectories on a population-based level and to illustrate which of different
personal and clinical characteristics are at highest risk of experiencing at
different time points following their cancer diagnosis.

### Author contributions

CEW, PCL and MJC contributed to the conception of the work. CEW did the data
collection. CEW and MJC implemented the methods, conducted the data analysis
and interpreted the findings. CEW, PWD and MJC drafted the article. All
authors made critical revision of the article and approved the version to be
published.

### Financial disclosure

None reported.

### Conflict of interest

The authors declare no potential conflict of interests.

## References

  * [1] Aziz N. M., Rowland J. H.. Trends and advances in cancer survivorship research: Challenge and opportunity. Seminars in Radiation Oncology. 2003;13:248–266.
  * [2] Floyd J. D., Nguyen D. T., Lobins R. L., Bashir Q., Doll D. C., Perry M. C.. Cardiotoxicity of cancer therapy. Journal of Clinical Oncology. 2005;23(30):7685–96.
  * [3] Lee M. S., Finch W., Mahmud E.. Cardiovascular complications of radiotherapy. American Journal of Cardiology. 2013;112(10):1688–96.
  * [4] Nimwegen F. A., Schaapveld M., Janus C. P., et al. Cardiovascular disease after Hodgkin lymphoma treatment: 40-year disease risk. JAMA Intern Med. 2015;175(6):1007–17.
  * [5] Huszti E., Abrahamowicz M., Alioum A., Binquet C., Quantin C.. Relative survival multistate Markov model. Statistics in Medicine. 2012;31:269–286.
  * [6] Gilard-Pioc S., Abrahamowicz M., Mahboubi A., et al. Multi-state relative survival modelling of colorectal cancer progression and mortality. Cancer Epidemiology. 2015;39:447–455.
  * [7] Gillaizeau F., Dantan E., Giral M., Foucher Y.. A multistate additive relative survival semi-Markov model. Statistical Methods in Medical Research. 2017;26:1700–1711.
  * [8] Sarfati D., Blakely T., Pearce N.. Measuring cancer survival in populations: relative survival vs cancer-specific survival. International Journal of Epidemiology. 2010;39:598–610.
  * [9] Dickman P. W., Sloggett A., Hills M., Hakulinen T.. Regression models for relative survival. Statistics in Medicine. 2004;23:51–64.
  * [10] Dickman P. W., Adami H.-O.. Interpreting trends in cancer patient survival. Journal of Internal Medicine. 2006;260:103–117.
  * [11] Weibull C. E., Björkholm M., Glimelius I., et al. Temporal trends in treatment-related incidence of diseases of the circulatory system among Hodgkin lymphoma patients.. International journal of cancer. 2019;.
  * [12] Lambert P. C., Dickman P. W., Nelson C. P., Royston P.. Estimating the crude probability of death due to cancer and other causes using relative survival models. Statistics in Medicine. 2010;29(7-8):885–895.
  * [13] Eloranta S., Lambert P. C., Andersson T. M. L., et al. Partitioning of excess mortality in population-based cancer patient survival studies using flexible parametric survival models. BMC Medical Research Methodology. 2012;12:86.
  * [14] Eloranta S., Lambert P. C., Sjöberg J., Andersson T. M. L., Bjorkholm M., Dickman P. W.. Temporal trends in mortality from diseases of the circulatory system after treatment for Hodgkin lymphoma: a population-based cohort study in Sweden (1973 to 2006). Journal of Clinical Oncology. 2013;31(11):1435–41.
  * [15] Putter H., Fiocco M., Geskus R. B.. Tutorial in biostatistics: competing risks and multi-state models. Statistics in Medicine. 2007;26:2389–2430.
  * [16] Iacobelli S., Carstensen B.. Multiple time scales in multi-state models. Statistics in Medicine. 2013;32:5315–5327.
  * [17] Crowther M. J., Lambert P. C.. Parametric multistate survival models: flexible modelling allowing transition-specific distributions with application to estimating clinically useful measures of effect differences. Statistics in Medicine. 2017;36:4719–4742.
  * [18] Fiocco M., Putter H., Houwelingen H. C.. Reduced-rank proportional hazards regression and simulation-based prediction for multi-state models. Statistics in Medicine. 2008;27:4340–4358.
  * [19] Gill Richard D, Johansen Soren, others . A survey of product-integration with a view toward application in survival analysis. The annals of statistics. 1990;18(4):1501–1555.
  * [20] Oakes David. Multiple time scales in survival analysis. Lifetime Data Analysis. 1995;1(1):7–18.
  * [21] Lee Minjung, Gouskova Natalia A, Feuer Eric J, Fine Jason P. On the choice of time scales in competing risks predictions.. Biostatistics (Oxford, England). 2017;18:15–31.
  * [22] Wolkewitz Martin, Cooper Ben S, Palomar-Martinez Mercedes, et al. Multiple time scales in modeling the incidence of infections acquired in intensive care units. BMC medical research methodology. 2016;16(1):116.
  * [23] Nelson C. P., Lambert P. C., Squire I. B., Jones D. R.. Flexible parametric models for relative survival, with application in coronary heart disease. Statistics in Medicine. 2007;26:5486–5498.
  * [24] Brenner H., Hakulinen T.. Substantial overestimation of standard errors of relative survival rates of cancer patients. American Journal of Epidemiology. 2005;161:781–786.
  * [25] Royston P., Parmar M. K. B.. Flexible parametric proportional-hazards and proportional-odds models for censored survival data, with application to prognostic modelling and estimation of treatment effects. Statistics in Medicine. 2002;21:2175–2197.
  * [26] Sjöberg J., Halthur C., Kristinsson S. Y., et al. Progress in Hodgkin lymphoma: A population-based study on patients diagnosed in Sweden from 1973-2009. Blood. 2012;119(4):990–6.
  * [27] Glimelius I., Ekberg S., Jerkeman M., et al. Long-term survival in young and middle-aged Hodgkin lymphoma patients in Sweden 1992-2009 – trends in cure proportions by clinical characteristics. American Journal of Hematology. 2015;90(12):1128–34.
  * [28] Maringe Camille, Belot Aurélien, Rubio Francisco Javier, Rachet Bernard. Comparison of model-building strategies for excess hazard regression models in the context of cancer epidemiology. BMC medical research methodology. 2019;19(1):210.
  * [29] Bower H., Andersson T. M-L, Crowther M. J., Dickman P. W., Lambe M., Lambert P. C.. Adjusting Expected Mortality Rates Using Information From a Control Population: An Example Using Socioeconomic Status.. American journal of epidemiology. 2018;187:828–836.
  * [30] Aalen O. O.. A linear regression model for the analysis of life times.. Statistics in medicine. 1989;8:907–925.
  * [31] Bluhmki Tobias, Putter Hein, Allignol Arthur, Beyersmann Jan, consortium COMBACTE-MAGNET. Bootstrapping complex time-to-event data without individual patient data, with a view toward time-dependent exposures.. Statistics in medicine. 2019;38:3747–3763.
  * [32] Allignol Arthur, Schumacher Martin, Wanner Christoph, Drechsler Christiane, Beyersmann Jan. Understanding competing risks: a simulation point of view. BMC medical research methodology. 2011;11(1):86.
  * [33] Hjort Nils Lid. On Inference in Parametric Survival Data Models. International Statistical Review / Revue Internationale de Statistique. 1992;60(3):355–387.

The methods developed in this paper have been implemented in the merlin and
multistate Stata packages, which can be installed using:

  * net install merlin, from(https://www.mjcrowther.co.uk/code/merlin)

  * net install multistate, from(https://www.mjcrowther.co.uk/code/multistate)

Further details and examples, on each package, can be found at
www.mjcrowther.co.uk/software/merlin and
www.mjcrowther.co.uk/software/multistate.

## Example Stata code

Here we detail example syntax that implements the proposed multi-state
modelling framework in Stata. Model estimates to obtain the predictions are
available upon request.

  *

        // Load analysis dataset
    use analysisdata, clear

    //===========================================================================================================//
    // Transition 1: Expected model

    // knots for age spline (log scale)
    local aknots 2.8904 4.0604 4.2195 4.3175 4.3944 4.5951
    // knots for year spline
    local yknots 1985 1991 1998 2003 2009 2014

    // estimate the spline-based Poisson model and store model object
    exptorcs                                              /// function name
                       female,                            /// baseline covariates in expected model
                       expdata(popinc)                    /// dataset with expected rates
                       event(_d) exposure(_y)             /// variables in exprates() file for poisson model
                       year(_year,                        /// year of diagnosis
                                 knots(`yknots´)               /// knot specification
                                 noorthog                      /// do not orthogonalise the spline terms
                       age(_age,                          /// age at diagnosis
                                 knots(`aknots´)               /// knot specification
                                 log                           /// create splines of log age
                                 noorthog)                     // do not orthogonalise the spline terms

    //store the estimated model
    estimates store expmod

    //===========================================================================================================//
    // We now reshape the analysis data into long format, with one row for each transition at which
    // each patient is at risk for

    // Using msset, which will assume a single starting state, and we enter variables
    // which contain event or censoring times for DCS or death. By default, an upper
    // triangular transition matrix will be assumed, i.e. the required illness-death model
    msset, id(lopnr) states(dcs dead) times(dcs_time dead_time)

    // now declare our dataset to be survival data
    stset _stop, enter(_start) failure(_status==1) scale(365.24)

    //===========================================================================================================//
    // Transition 2: Excess model

    merlin (_t female rcs(_year, knots(`yknots´))                            ///
                      rcs(_age, knots(`aknots´) log)                         ///
                      if _trans==1,                                          ///
                      family(rp, failure(_d) ltruncated(_t0)                 ///
                        df(5) bhazard(exp_dcsrate) noorthog))
    estimates store excessmod

    //===========================================================================================================//
    // Transition 3: Cancer diagnosis to death (no DCS)

    merlin (_t female rcs(_year, knots(`yknots´))                            ///
                      rcs(_age, knots(`aknots´) log)                         ///
                      if _trans==2,                                          ///
                      family(rp, failure(_d) ltruncated(_t0)                 ///
                        df(5) noorthog))
    estimates store deadmod

    //===========================================================================================================//
    // Transition 4: DCS to death

    // We must first generate the clock-reset timescale
    gen t2 = _stop - _start

    // and re-stset
    stset t2, failure(_status==1) scale(365.24)

    // estimate the model
    merlin (_t female rcs(_year, knots(`yknots´))                            ///
                      rcs(_age, knots(`aknots´) log)                         ///
                      if _trans==3,                                          ///
                      family(rp, df(4) noorthog))
    estimates store dcsdeadmod

    //===========================================================================================================//
    // Predictions

    // Define transition matrix of the proposed extended illness-death model with the expected and excess
    // components
    mat tmat =  (.,1,2,3,.\.,.,.,.,4\.,.,.,.,5\.,.,.,.,.\.,.,.,.,.)

    // Define age and year groups that we want predictions for
    local age 30
    local year 1995

    // Generate a time variable at which to obtain predictions (this will be on the main time since cancer
    // diagnosis timescale)
    range tvar 0 15 1000

    // Calculate transition probabilities for a male aged 30 at diagnosis, diagnosed in 1995
    predictms,   transmatrix(tmat)                 ///
                 models(expmod excessmod deadmod   /// pass model objects for transitions 1 to 5
                   dcsdeadmod dcsdeadmod)          ///
                 probability                       /// predict transition probabilities
                 timevar(tvar)                     /// time points to predict at
                 at1(female 0                      /// specify the covariate pattern to predict at
                   _year `year´                    /// specify y, year of diagnosis
                   _age `age´)                     /// specify a, age of diagnosis
                 tsreset(4 5)                      /// transitions 4 and 5 are on the clock-reset timescale
                 ci                                 // calculate confidence intervals

    // Calculate differences in transition probabilities between females and males
    // aged 30 at diagnosis, diagnosed in 1995
    predictms,   transmatrix(tmat)                 ///
                 models(expmod excessmod deadmod   /// pass model objects for transitions 1 to 5
                   dcsdeadmod dcsdeadmod)          ///
                 probability                       /// predict transition probabilities
                 timevar(tvar)                     /// time points to predict at
                 at1(female 0                      /// specify the 1st covariate pattern to predict at
                    _year `year´                    /// specify y, year of diagnosis
                   _age `age1´)                    /// specify a, age of diagnosis
                 at2(female 1                      /// specify the 2nd covariate pattern to predict at
                   _year `year´                    /// specify y, year of diagnosis
                   _age `age´)                     /// specify a, age of diagnosis
                 tsreset(4 5)                      /// transitions 4 and 5 are on the clock-reset timescale
                 difference                        /// calculate the difference in trans. probs. between at2() - at1()
                 ci                                 // calculate confidence intervals

## Supplementary figures

SUPPLEMENTARY FIGURE 8: Probability of being in each state for male patients
diagnosed in 1995, aged 30 years, as a function of time since diagnosis as
predicted from 18 different models with varying degrees of freedom. State
definitions: 1 - HL, 2 - Excess DCS, 3 - Expected DCS, 4 - Dead before DCS, 5
- Dead after DCS. SUPPLEMENTARY FIGURE 9: Probability of being in each state
for male patients diagnosed in 1995, aged 60 years, as a function of time
since diagnosis as predicted from 18 different models with varying degrees of
freedom. State definitions: 1 - HL, 2 - Excess DCS, 3 - Expected DCS, 4 - Dead
before DCS, 5 - Dead after DCS.


Using The Newton-Raphson method with Automatic Differentiation to numerically
solve Implied Volatility of stock option through Binomial Model Wanchaloem
Wunkaew 1,2 , Yuqing Liu 1 and Kirill V. Golubnichiy 1 1 Department of
Mathematics, University of Washington, Seattle, WA 98195 2 Department of
Applied Mathematics, University of Washington, Seattle, WA 98195 Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the paper written by Klibanov et al [3], it proposes a novel method to
calculate implied volatility of a European stock options as a solution to ill-
posed inverse problem for the Black-Scholes equation. In addition, it proposes
a trading strategy based on the difference between implied volatility of the
option and the volatility of the underlying stock. In addition to the Black-
Scholes equation, Binomial model is another method used to price European
options. And, the implied volatility can be also calculated through this
model. In this paper, we apply the Newton-Raphson method together with
Automatic Differention to numerically approximate the implied volatility of an
arbitrary stock option through this model. We provide an explanation of the
mathematical model and methods, the methodology, and the results from our test
using the stimulated data from the Geometric Brownian Motion Model and the
Binomial Model itself, and the data from the US market data from 2018 to 2021.

Keywords: Binomial Model; Implied Volatility; Option Pricing; Black-Scholes
Equation; Automatic Differentiation; Newton-Raphson Method.

##  1 Introduction

An option is a contract that allows the holder the right to buy or sell an
underlying asset or financial instrument at a specific strike price on or
before a specific date. An option to buy is known as a ”call option,” and an
option to sell is known as a ”put option.” Options are divided into two types
based on the type of time period of implementation. An American type of option
is an option that can be exercised at or before maturity, whereas a European
option is an option that can only be exercised at the maturity date. In this
paper, we used the market data for a total of 30,000 European call stock
options.

The Black-Scholes model, also known as the Black-Scholes-Merton model, is one
of the most widely used models in modern financial theory when estimating the
theoretical value of options. It was developed in 1973 by Fischer Black,
Robert Merton, and Myron Scholes and was initially introduced in [1]. In the
model, it uses the current stock prices, expected dividends, the option’s
strike price, expected interest rates, time to expiration, and expected
volatility to calculate the theoretical value of an option contract. This
equation can be written as a parabolic partial differential equation which can
be explicitly solved. Despite its simplicity, solving implied volatility via
the equation is difficult due to its ill-posedness characteristic [3].

Another commonly-used model for option valuation is the binomial option
pricing model [2]. It was first proposed by William Sharpe in the 1978 and
then formalized by Cox, Ross, and Rubinstein in 1979. It traces the evolution
of an option’s key underlying variables in discrete time by means of a
binomial tree. In this paper, we will discuss the method of using the binomial
option pricing model to calculate the implied volatility. We will explain the
methodology behind the idea as well as some numerical results from our
experiments.

In the paper [3], it proposed a new mathematical method is accurately
calculate the implied volatility of a European stock option. In addition, it
derives a trading strategy based on the difference between the implied
volatility of an option and the volatility of its underlying stock.
Specifically, let $\sigma$ be the volatility of a stock, $\hat{\sigma}$ be the
volatility of the corresponding option to the stock. Let $S$ be the price of
the stock, $f(S,K)$ be the payoff or value function of the option at the
maturity time T, and $u(S,t)$, with t, the time until the maturity will occur,
be the price of the option.

By assumptions from the paper, the following equation is satisfied [3]:

$\displaystyle\frac{\partial u(S,t)}{\partial t}$
$\displaystyle=\frac{\hat{\sigma}^{2}}{2}S^{2}\frac{\partial^{2}u(S,t)}{\partial
S^{2}},\qquad S>0$ $\displaystyle u(S,0)$ $\displaystyle=f(S)$

Let K be the strike price, then the payoff function is $f(S)=\max(S-K,0)$, and
$u(S,t)$ is given by the Black-Scholes formula [1]:

$\displaystyle u(S,t)=S\Phi(\Theta_{+}(S,t))-e^{-rt}K\Phi(\Theta_{-}(S,t))$

where the risk-free rate $r$ is assumed to be $0$. The stochastic equation of
the Geometric Brownian motion for the stock price $S_{t}$ with the volatility
$\sigma$ and zero drift is expressed by

$\displaystyle dS_{t}=\sigma S_{t}dW_{t}$

Since the expectation of the Wiener process, dW, is 0 and that by substituting
the above equations in the Itô formula, we get the expected value of the
increment of the option price on an infinitesimal time interval is

$\displaystyle\frac{(\sigma^{2}-\hat{\sigma}^{2})}{2}S^{2}\frac{\partial^{2}u(S,t)}{\partial^{2}S}dt$

Since $\frac{\partial^{2}u(S,t)}{\partial^{2}S}$ is the Greek $\Gamma$ which
is proved to be non-negative, and $S^{2}$ is greater than $0$; thus, the
expected direction of an option price depends on $\sigma-\hat{sigma}$; that
is, under the assumptions, we expect the option price in the market to
increase if $\sigma-\hat{sigma}>0$, the volatility of the underlying asset is
greater than the implied volatility of the option. After deriving the
equation, the paper also examines the validity of the mathematical model and
generates a possible winning strategy from the model. After testing the
hypothesis, it gets a precision of 55.57%.

The result from [3] strengthens the role and the importance of the implied
volatility of options in financial mathematics. However, the Black-Scholes
equation, on which the paper relies on, is not the only single method used to
price an option. Moreover, this method works only for European options.

In this paper, we firstly present the algorithm of Binomial model and explain
why it can be used to find the implied volatility of an option’s price by
showing the convergence from the Black-Scholes equation to the Binomial model.
We discuss the method of calculating implied volatility inversely via binomial
model by explaning the conditions and parameters the model is required.

Then, we test our model by employing the binomial model to calculate the
implied volatility of an option price generated from Geometric Brownian motion
model with $\mu=0$, $\sigma=0.2$, and compare the result from the Binomial
Model to the result from the Black-Scholes Model, and find out the relative
difference between the precision of two models.

##  2 Mathematical Model

###  2.1 Binomial Model

The Binomial Model is based on the main assumption that the price of the
underlying asset, which is the stock in this case, can change in only two
possible ways: either going up or going down, in the next time step. These
stock prices can be estimated using Geometric Brownian motion:

$\displaystyle dS_{t}$ $\displaystyle=\mu S_{t}dt+\sigma S_{t}dW_{t}$

where $S_{t}$ is a stock price at time $t$, $\mu$ is a drift, $\sigma$ is an
volatility of the stock, and $W_{t}$ is a Wiener process.
Note that in the next step $\delta t$, the the Geometric Brownian motion can
be approximated into two possible values [6]:

$\displaystyle S_{t+\delta t}$ $\displaystyle=S_{t}e^{\sigma\sqrt{\delta
t}}\quad\text{ or }\quad S_{t+\delta t}=S_{t}e^{-\sigma\sqrt{\delta t}}$

Assume that the next time step, at time $T$, is the maturity of a call option
with strike price $K$, the value of the function in that step can be
determined by the function:

$\displaystyle V(T)$ $\displaystyle=\max(0,S(T)-K)$

where $S(T)$ is the price of the stock at time $T$ (the maturity of the call
option).

Since we assume that the price of the stock can go in only two directions in
the next step, we can separately calculate the value of the call option as
follows:

$\displaystyle C^{+}$ $\displaystyle=\max(0,S(0)e^{\sigma\sqrt{\delta
t}}-K)=\max(0,S_{0}e^{\sigma\sqrt{T}}-K)$ $\displaystyle C^{-}$
$\displaystyle=\max(0,S(0)e^{-\sigma\sqrt{\delta
t}}-K)=\max(0,S_{0}e^{-\sigma\sqrt{T}}-K)$

Then, A single-step Binomial Model is formulated as shown below:

$\displaystyle C=e^{-r{\delta}t}\left((\frac{1}{2}+\frac{r\sqrt{\delta
t}}{2\sigma})C^{+}+(\frac{1}{2}-\frac{r\sqrt{\delta t}}{2\sigma})C^{-}\right)$

where $C$ is the value of option; $r$ is a risk-free rate; $t$ is a time-step
(in year); $\sigma$ is a volatility of the underlying asset.

This model can be extended into multiple steps, which is called the Lattice
Model. In this case, the values of options in any possible branch of the tree
can be determined by the value functions as shown above. Yet, in the middle
steps, the values of options in the steps after them are required in order to
calculate the option at that middle step. This can be done by solving the
option values in each step backwards. In the last step, the value of the
option at the current time is determined.

###  2.2 Transformation/convergence from Black-Scholes to Binomial Model

The fundamental difference between the Black-Scholes equation and the binomial
model is that Black-Scholes equation considers the continuous time, whereas
the binomial model considers time as a discrete variable.

One way of deriving Black-Scholes equation is letting the number of time steps
in the binomial tree to approach infinity.

Suppose a binomial model with n time steps is used to value a European call
option, then each step has the length T/n. It would have i upward movements
and n-i downward movements. According to the binomial model, the final stock
price is $S(0)u^{i}d^{n-i}$, where S(0) is the original price, u is the up
factor, and d is the down factor.

In order to earn a profit from the option, we make a move when
$\max(0,S(0)u^{i}d^{n-i}-K)>0$. In other words, $S(0)u^{i}d^{n-i}-K>0$, and
$S(0)u^{i}d^{n-i}>K$.

By substituting the equation for u and d, we get:

$\displaystyle\ln(\frac{S(0)}{K})>n\sigma\sqrt{\delta t}-2i\sigma\sqrt{\delta
t}$

or

$\displaystyle i>\frac{n}{2}-\frac{\ln(\frac{S(0}{K})}{2\sigma\sqrt{\delta
t}}$

Then, the equation of C can be written as:

$\displaystyle C=e^{-rT}(S(0)A_{1}-KA_{2})$

where

$\displaystyle
A_{1}=\sum_{i>\frac{n}{2}-(ln(\frac{S_{0}}{K})2\sigma\sqrt{\delta{t}}}\frac{n!}{(n-1)!i!}p^{i}(1-p)^{n-1}u^{i}d^{n-1}$

$\displaystyle
A_{2}=\sum_{i>\frac{n}{2}-(ln(\frac{S_{0}}{K})2\sigma\sqrt{\delta{t}}}\frac{n!}{(n-1)!i!}p^{i}(1-p)^{n-1}$

Since the binomial distribution approaches a normal distribution as the number
of trials approaches infinity, and the mean is np and standard deviation is
np(1-p), with p the probability of success (in this case, the probability of
going up), n the number of trials.

Then

$\displaystyle
A_{2}=N(\frac{np-(\frac{n}{2}-(ln(\frac{S(0)}{K}){2\sigma\sqrt{\delta{t}}})}{\sqrt{np(1-p)}})=N(\frac{ln(\frac{S(0)}{K})+(r-\frac{\sigma^{2}}{2})\delta{t}}{\sigma\sqrt{\delta{t}}})$

as n goes to infinity. Similarly for A1, we get

$\displaystyle
A_{1}=e^{rt}N(\frac{ln(\frac{S(0)}{K})+(r+\frac{\sigma^{2}}{2})\delta{t}}{\sigma\sqrt{\delta{t}}})$

Combining $A_{1}$, $A_{2}$, and C, we get

$\displaystyle
d_{1}=\frac{\ln(\frac{S(0)}{K})+(r+\frac{\sigma^{2}}{2}\delta{t})}{\sigma\sqrt{\delta{t}}}$

$\displaystyle
d_{2}=\frac{\ln(\frac{S(0)}{K})-(r+\frac{\sigma^{2}}{2}\delta{t})}{\sigma\sqrt{\delta{t}}}=d_{1}-\sigma\sqrt{T}$

which is the Black-Scholes equation for finding a European call option.

###  2.3 Numerical Methods

####  2.3.1 Newton-Raphson method

The Newton-Raphson method is widely used to numerically solve roots of
differentiable functions [7]. The algorithm can be written as the formula
shown below:

$\displaystyle x_{n+1}$
$\displaystyle=x_{n}-\frac{f(x_{n})}{f^{\prime}(x_{n})}$

where $x_{n}$ for $n\in\mathbb{N}$ is an approximated root for $n^{\text{th}}$
step, and $x_{0}$ is an initial guess, and $f$ is a function. Under stable
conditions, $x_{n}$ will approach to one root of $f$.

####  2.3.2 Automatic Differentiation

As shown in the previous subsection, the Newton-Raphson method requires a
derivative of $f$ in order to calculate a root for $f$. There are plenty of
methods to approximate the derivative. One of the intuitive methods is to
directly calculate it using theorems from calculus, so-called symbolic
differentiation. However, because some formulas, such as product rules and
chain rules, greatly expand functions, this method may result in an overly
complex expression. This symbolic method may be applicable for single-step
binomial models, but it is exhaustive for multi-step lattice models. Another
numerical method which is widely used is the finite difference method:

$\displaystyle f^{\prime}(x)$ $\displaystyle\approx\frac{f(x+\Delta
x)-f(x)}{\Delta x}$

for some small $\Delta x>0$. The finite element method is also introduced to
calculate the derivative for the Newton-Raphson method, and the modified
method is called the Secant Method.

In this project, we apply a widely-used method for differentiation, so-called
Automatic Differentiation (AD). AD calculates the derivative of a function
immediately after each component of the function has been computed. Also, the
derivative of that part of the function is stored together with the result of
that part being stored[5].

Automatic Differentiation is largely used in machine learning algorithms such
as back propagation. With the structure of the Binomial/Lattice Model, the
option price is calculated as a multiple of steps of calculation through a
tree. We propose to use Automatic Differentiation for calculating the
derivative of the model function.

##  3 Methodology

In this project, we approximate the implied volatility of a European stock
option by using the Binomial/Lattice model.

###  3.1 Calculating implied volatility inversely via binomial model

As mentioned in the previous section, in order to determine the current value
of an option, some parameters such as the current stock price and risk-free
rate are required.These parameters can be obtained from the exchange market.
Moreover, in the same place, the price of the option is determined.

Now, we let the known parameters be constants. Then, we can simply write a
formula as:

$\displaystyle C$ $\displaystyle=f(\sigma)$

Then, the $\sigma$ satisfying the equation above is the implied volatility of
a corresponding stock option. It is important to note that in this paper, the
10-layer lattice model is applied. The number of the layers is determined by
computational limitations. The risk-free rate is based on the 10-year treasury
yield for the real-market data. And we assume the rate to be $0$for the
generated dataset. Other parameters are obtained from the real market. C is
known in the exchange market. Thus, $\sigma$ can be determined by finding the
root of the equation:

$\displaystyle f(\sigma)-C$ $\displaystyle=0$

We denote the root of the equation above by $\hat{\sigma}$.
To finding a root of the equation, we apply Newton-Raphson’s method with
initial value, $x_{0}=0.2$, with the tolerant rate $=10^{-5}$. That is, the
iteration would stop at $n$ when $|f(x_{n})-f(x_{n-1})|<10^{-5}$. The maximum
number of iterations is set to 100. The derivative of $f$ is calculated using
Automatic Differentiation in each iteration.

###  3.2 Testing the numerical method

Since the model used in this project is combined with multiple step
calculations, it is difficult for the numerical method mentioned above to
accurately approximate the implied volatility. Therefore, to assess the
efficiency of the model, we tested the algorithms with simulated data.

Similar to [3], the underlying stock is simulated using a Geometric Brownian
motion model with $\mu=0$ and $\sigma=0.2$. And, the price of options are
calculated via the same binomial model, with 90 days of maturity and with 33
different volatilities. After the simulated option price is obtained, we
inversely solve for the implied volatility with the method mentioned above.

Figure 1: Simulated stock price based on Geometric Brownian motion

Figure 2: An example of simulated option prices generated from the stock
prices above and Binomial model with four different volatility , strike price
of 150, and fixed maturity at day 180.

###  3.3 Testing with the real world data

After testing the method with synthetic data, we assess the efficiency of the
model with real-world data. In this project, we apply the algorithm to 38270
European options traded in the market from 2018 to 2021. The last prices of an
option and its underlying asset are used in this evaluation.

##  4 Numerical Finding

The scatter plots show the implied volatility from the binomial model with
their corresponding actual implied volatility. The red diagonal line plot
shows the position when a pair of implied volatility are the same. Note that
the algorithm could possibly return negative implied volatility, which is not
shown in these plots due to its absurdity.

According to the plots, the algorithm frequently underestimates the actual
implied volatility. To address the issues of underestimation and negative
results, we show the histograms together with the scatter plot for each
dataset.

###  4.1 Generated Data

Figure 3: (left) The scatter plot of implied volatilities calculated from the
binomial model and their corresponding actual values from the dataset
(Generated Data); (right) The implied volatilities distributions (Generated
Data).

###  4.2 Stock options selected from American markets from 2017-2018

Figure 4: (left) The scatter plot of implied volatilities calculated from the
binomial model and their corresponding actual values from the dataset (The
2017-2018 US market); (right) The implied volatilities distributions (The
2017-2018 US market).

###  4.3 Stock options selected from American markets from 2020-2021

Figure 5: (left) The scatter plot of implied volatilities calculated from the
binomial model and their corresponding actual values from the dataset (The
2020-2021 US market); (right) The implied volatilities distributions (The
2020-2021 US market).

##  5 Discussion

The three results show this method greatly underestimates the actual implied
volatility. Among the plots, the result from testing the method on generated
data best shows the efficiency of the method. Yet, testing with the real-world
dataset reveals some flaws in its efficiency, particularly in the 2020-2021
dataset. One of the important causes of this is that the collection of actual
implied volatilities attained in the dataset was inversely solved from Black-
Scholes Models. This results in high error rates in real-world data but much
less in generated data. Comparing the results from different periods of real
market data, we found that the results from 2020 to 2021 attained a higher
error than those of the 2017-2018 market. This might be from the effect of the
pandemic on the option market.
Among the experiments, the results from applying the proposed method to the
generated data are the most acceptable. The absolute difference between the
actual value and the result from the method is less than $0.01$ for each
simulated option price. However, most of the results underestimate the actual
implied volatility, which is shown on the corresponding histogram. Despite the
fact that it is not clear how this occurs, the distribution of the difference
seemingly allows us to adjust the results from the algorithm by adding the
constant average error from the dataset to gain more accurate results. Yet,
the adjustment needs to be deeply studied in the future.
The result from applying the method to the real stock option market from 2017
to 2018 shows more error compared to that on the generated data. In addition,
the method does not wholly underestimate the actual values; it sometimes
overestimates the values. However, the distribution plot still implies that
the method averagely underestimates the values.

Another interesting point is that the numerical method gives negative implied
volatility, which is absurd. We provide three explanations for these issues:

  1. 1.

European stock options are priced based on the Black-Scholes equation, not the
Binomial Model. Despite the fact that the Binomial Model converges to the
model in continuous time, the implied volatility inversely calculated from the
different models may not necessarily be the same. Presumably, there is no
positive or real root, or implied volatility, for the Binomial Model.

  2. 2.

The Newton-Raphson method gives only one root. Even if we assign the initial
value,$x_{0}$, to be close to the average option price, under some unstable
conditions, the Newton-Raphson method may give other possible roots, which are
not necessarily negative, of the function. This issue can be fixed by testing
the method on the different initial values and selecting the most appropriate
result.

  3. 3.

Unobservable real-world factors. Unlike the generated dataset, where
everything is under control, the real stock market is based on uncertainty.

The results from applying the method to the real option market dataset from
2020 to 2021 show the largest error. As shown in the scatter plot, the method
greatly underestimates the implied volatility. The explanations from the
previous paragraphs can account for most of the results on this dataset.
Importantly, the uncertainty of the market during the pandemic may influence
the results.

##  6 Conclusion

In this paper, we provide a method to numerically calculate the implied
volatility of European stock options. The method consists of two parts: the
Binomial Model and a numerical method. The Binomial model is used to price
European stock options. Given a data about options in the exchange market,
including the price of the options, we can calculate the implied volatility by
inversely solving the model to get a volatility variable pertaining to the
given option price. Newton-Raphson method is used to solve for such
volatility. Yet, the numerical method requires the derivative of the binomial
model, which is hard and tedious to analytically or numerically determine. We
propose to use Automatic Differentiation, which immediately can calculate the
derivative of function after the function parameters are passed. Finally, we
test the proposed method on three datasets: the generated data, the real-
market dataset from 2017 to 2018, and that from 2020 to 2021. The results show
that the method works well for the generated data, but its efficiency deceases
once tested with the real-market data; yet, the result from 2017-2018 dataset
is acceptable. We additionally provide the explanations of the issues found in
testing the method with each dataset.
Practically, we can use the calculated implied volatility to make an
interference in option trading decision as provided in. In addition, this
method can be modified in order to numerically calculate some parameters for
other asset pricing models. For example, this method can be applied to
calculate the implied volatility for an American option which does not follow
the Black-Scholes equation; yet, this kind of option can be calculated using
modified binomial model.

## References

  * [1] Black, F., & Scholes, M. (2019). The pricing of options and corporate liabilities. In World Scientific Reference on Contingent Claims Analysis in Corporate Finance: Volume 1: Foundations of CCA and Equity Valuation (pp. 3-21).
  * [2] Cox, J. C., Ross, S. A., & Rubinstein, M. (1979). Option pricing: A simplified approach. Journal of financial Economics, 7(3), 229-263.
  * [3] Klibanov, M. V., Shananin, A. A., Golubnichiy, K. V., & Kravchenko, S. M. (2022). Forecasting Stock Options Prices via the Solution of an Ill-Posed Problem for the Black-Scholes Equation. arXiv preprint arXiv:2202.07174.
  * [4] Kim, Y. S., Stoyanov, S., Rachev, S., & Fabozzi, F. (2016). Multi-purpose binomial model: Fitting all moments to the underlying geometric Brownian motion. Economics Letters, 145, 225-229.
  * [5] Rall, L. B., & Corliss, G. F. (1996). An introduction to automatic differentiation. Computational Differentiation: Techniques, Applications, and Tools, 89, 1-18.
  * [6] Ross, S. M. (2011). An elementary introduction to mathematical finance. Cambridge University Press.
  * [7] Corless, R. M., & Fillion, N. (2013). A graduate introduction to numerical methods. AMC, 10, 12.


SNSN-323-63

Direct CPV in two-body and multi-body charm decays at LHCb

Evelina Gersabeck on behalf of the LHCb collaboration

Ruprecht-Karls-Universitaet Heidelberg
Physikalisches Institut
Im Neuenheimer Feld 226
69120 Heidelberg, Germany

> The Standard Model predicts CP asymmetries in charm decays of $O(10^{-3})$
> and the observation of significantly larger CP violation could indicate non-
> Standard Model physics effects. During 2011 and 2012, the LHCb experiment
> collected a sample corresponding to $3/fb$ yielding the world s largest
> sample of decays of charmed hadrons. This allowed the CP violation in charm
> to be studied with unprecedented precision in many two- body and multibody
> decay modes. The most recent LHCb searches for direct CP violation are
> presented in these proceedings.

> PRESENTED AT
>
>
>
>
> Presented at the 8th International Workshop on the CKM Unitarity Triangle
> (CKM 2014), Vienna, Austria,
> September 8-12, 2014

##  1 Introduction

The excellent performance of the LHC and the LHCb experiment, along with large
production $c\overline{c}$ cross sections for $pp$ collisions at $\sqrt{s}$ of
7 and 8 TeV has enabled unprecedentedly large samples of charm decays to be
recorded during 2011 and 2012, corresponding to $3/fb$ of integrated
luminosity. These large samples allow the study of CP violation (CPV) effects
at a precision not achieved before in charm decays. The data was taken with a
regular swap of the polarity of the spectrometer dipole magnet which can
compensate for the left-right detector asymmetries to a first order. Both
charm decays, promptly produced in the primary pp interaction, and coming from
a parent beauty hadron are exploited at LHCb; this is indicated for each of
the presented analyses.

##  2 Time-integrated CP asymmetry in ${{\cal D}^{0}}\rightarrow h^{+}h^{-}$
from semileptonic decays

A search for a time-integrated CP asymmetry in ${{\cal D}^{0}}\rightarrow
h^{+}h^{-}$ decays is performed using the full dataset corresponding to
$3/fb$. The flavour of the initial ${{\cal D}^{0}}$ state is tagged by the
charge of the muon in the semileptonic $B\rightarrow{{\cal
D}^{0}}\mu^{-}\nu_{\mu}X$ decays.

The raw measured asymmetry for tagged ${{\cal D}^{0}}$ messns to a final state
$f$ is given by:

$A_{raw}(f)=\frac{N(B\rightarrow{{\cal D}^{0}}\mu^{-}X)-N(B\rightarrow{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}\mu^{+}X)}{N(B\rightarrow{{\cal
D}^{0}}\mu^{-}X)+N(B\rightarrow{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}\mu^{+}X)},$ (1)

where $N$ indicates the number of reconstructed events of a given decay after
background subtraction, and $X$ refers to the undetected final state particles
from the semileptonic B decay. The raw asymmetry is a sum of the physical CP
asymmetry, ($A_{CP}(f)$), the production asymmetry ($A_{P}(B)$) and detection
asymmetry ($A_{D}(\mu)$) :

$A_{raw}(f)=A_{CP}(f)+A_{P}(B)+A_{D}(\mu).$ (2)

As the quantity of interest is $A_{CP}(f)$, the main experimental challenge is
to separate it from the nuisance asymmetries. An experimentally more robust
variable, $\Delta A_{CP}$ can be constructed by taking the difference of the
raw asymmetries measured in ${{\cal D}^{0}}\rightarrow K^{+}K^{-}$ and ${{\cal
D}^{0}}\rightarrow\pi^{+}\pi^{-}$ decays:

$\Delta A_{CP}=A_{raw}(KK)-A_{raw}(\pi\pi),$ (3)

and thus cancelling the production and the muon detection asymmetries to a
first order. Alternatively, for extracting $A_{CP}(KK)$, the detection and
production asymmetries can be measured using Cabibbo-favoured (CF)
$B\rightarrow{{\cal D}^{0}}(\rightarrow K^{-}\pi^{+})\mu^{-}X$ decays where no
CPV is expected. An additional detection asymmetry, $A_{D}(K\pi)$, arises due
to the different interaction rates of the charged $K$ with the matter. To
remove this asymmetry, the control channels $D^{+}\rightarrow
K^{-}\pi^{+}\pi^{+}$ and $D^{+}\rightarrow\overline{K^{0}}\pi^{+}$ are used.
In the $D^{+}\rightarrow\overline{K^{0}}\pi^{+}$ decays, the detection
asymmetry of $K^{0}$ arising due to the combined effect of CPV in mixing in
the neutral kaon system and the different interaction rates of $K^{0}$ and
$\overline{K^{0}}$ in the detector material, $A_{CP/int}$, is estimated from
simulation and subtracted from the raw asymmetry. Once $\Delta A_{CP}$ and
$A_{CP}(KK)$ are measured, the individual asymmetry $A_{CP}(\pi\pi)=\Delta
A_{CP}-A_{CP}(KK)$ can be computed.

In total, $\sim 2.1\mathrm{\times}10^{6}~{}{{\cal D}^{0}}\rightarrow KK$ and
$\sim 0.7\times 10^{6}~{}{{\cal D}^{0}}\rightarrow\pi\pi$ decays are
reconstructed. The analysis is done separately for the 2011 and the 2012 data,
and for the two magnet polarities and consistent results are obtained. As the
production and detection asymmetries depend on the kinematic distributions,
weights are assigned to the candidates such that the kinematic distributions
are equalised. The results for the asymmetries [1]

$\Delta A_{CP}=(+0.14\pm 0.16(stat)\pm 0.08(syst))\%,$ (4)
$A_{CP}(KK)=(-0.06\pm 0.15(stat)\pm 0.10(syst))\%,$ (5)

and, with correlation rho=0.28,

$A_{CP}(\pi\pi)=(-0.20\pm 0.19(stat)\pm 0.10(syst))\%$ (6)

are compatible with CP conservation. The $A_{CP}(hh)$ asymmetries are the most
precise measurement of individual asymmetries up to date. The precision of
$\Delta A_{CP}$ is comparable to the preliminary result for $\Delta A_{CP}$
measured using prompt ${{\cal D}^{0}}$ decays reconstructed in 1/fb of
integrated luminosity [2] .

##  3 CP asymmetries in $D^{+}_{(s)}\rightarrow K^{0}_{S}h^{+}$ decays

Following a similar analysis strategy, the raw asymmetries in the prompt
$D^{+}_{(s)}\rightarrow K^{0}_{S}h^{+}$ decays

$A_{raw}(K^{0}_{S}h^{+})=A_{CP}(K^{0}_{S}h^{+})+A_{P}(D^{+}_{(s)})+A_{D}(h^{+})+A_{CP/int}(K^{0}/\overline{K^{0}}),$
(7)

where $h^{+}=K^{+}$ or $\pi^{+}$, are used to extract the quantities of
interest

$\displaystyle A_{CP}^{DD}$ $\displaystyle=$
$\displaystyle((A_{raw}(D^{+}_{s}\rightarrow
K^{0}_{S}\pi^{+})-(A_{raw}(D^{+}_{s}\rightarrow K^{0}_{S}K^{+}))$
$\displaystyle-$ $\displaystyle((A_{raw}(D^{+}\rightarrow
K^{0}_{S}\pi^{+})-(A_{raw}(D^{+}\rightarrow K^{0}_{S}K^{+}))$
$\displaystyle\approx$ $\displaystyle A_{CP}(D^{+}_{s}\rightarrow
K^{0}_{S}\pi^{+})+A_{CP}(D^{+}\rightarrow K^{0}_{S}K^{+}),$

and

$A_{CP}(D^{+}_{S}\rightarrow K^{0}_{S}\pi^{+})=A_{raw}(D^{+}_{s}\rightarrow
K^{0}_{S}\pi^{+})-A_{raw}(D^{+}_{s}\rightarrow\Phi\pi^{+}),$ (9)
$\displaystyle A_{CP}(D^{+}\rightarrow K^{0}_{S}K^{+})$ $\displaystyle=$
$\displaystyle(A_{raw}(D^{+}\rightarrow
K^{0}_{S}K^{+})-A_{raw}(D^{+}_{s}\rightarrow K^{0}_{S}K^{+}))$
$\displaystyle-$ $\displaystyle(A_{raw}(D^{+}\rightarrow
K^{0}_{S}\pi^{+})-A_{raw}(D^{+}_{S}\rightarrow\Phi\pi^{+})).$

The analysis is done using the full available data sample of $3/fb$. In total,
$\sim 4.8\times 10^{6}~{}D^{+}\rightarrow K_{S}^{0}\pi^{+}$, $\sim 0.12\times
10^{6}~{}D^{+}_{S}\rightarrow K_{S}^{0}\pi^{+}$, $\sim 1.0\times
10^{6}~{}D^{+}\rightarrow K_{S}^{0}K^{+}$, $\sim 1.5\times
10^{6}~{}D^{+}_{S}\rightarrow K_{S}^{0}K^{+}$, $\sim 7\times
10^{6}~{}D^{+}\rightarrow\Phi\pi^{+}$, and $\sim 13.\times
10^{6}~{}D^{+}_{S}\rightarrow\Phi\pi^{+}$ decays were reconstructed. The
analysis is done separately for 2011 and 2012 data, and for both magnet
polarities and consistent results are obtained.. The results

$A_{CP}(D^{+}\rightarrow K^{0}_{S}K^{+})=(+0.03\pm 0.17(stat)\pm
0.14(syst))\%,$ (11) $A_{CP}(D^{+}_{S}\rightarrow K^{0}_{S}\pi^{+})=(+0.38\pm
0.46(stat)\pm 0.17(syst))\%,$ (12) $A_{CP}(D^{+}\rightarrow
K^{0}_{S}K^{+})+A_{CP}(D^{+}_{S}\rightarrow K^{0}_{S}\pi^{+})=(0.41\pm
0.49(stat)\pm 0.26(syst))\%.$ (13)

show no indication of CPV [3]. These are the most precise measurements of
these quantities.

##  4 Search for CPV in ${{\cal D}^{0}}\rightarrow\pi^{-}\pi^{+}\pi^{0}$
decays with the energy test

The energy test [4] is an unbinned model-independent statistical method to
search for time- integrated CP violation in ${{\cal
D}^{0}}\rightarrow\pi^{-}\pi^{+}\pi^{0}$ decays. The method relies on the
comparison of two ${{\cal D}^{0}}$ and $\overline{{{\cal D}^{0}}}$ flavour
samples and is sensitive to CPV localised in the phase-space of the multi body
final state. The flavour of the prompt ${{\cal D}^{0}}$ is tagged by the
charge of the slow pion in the decay $D^{*}\rightarrow{{\cal D}^{0}}\pi_{s}$.
For the reconstruction of the ${{\cal D}^{0}}$ both merged and resolved
neutral pions are used. The previous most sensitive study of this decay has
been done by the BaBar collaboration [5]. At LHCb, the energy test is used to
assign a p-value for a non-zero CPV hypothesis [6]. In this method, a test
statistic $T$ is used to compare the average distances based on the metric
function $\psi$. It is defined as

$T=\sum_{i,j>i}^{n}\frac{\psi_{ij}}{n(n-1)}+\sum_{i,j>i}^{\overline{n}}\frac{\psi_{ij}}{\overline{n}(\overline{n}-1)}-\sum_{i,j}^{n,\overline{n}}\frac{\psi_{ij}}{n\overline{n}},$
(14)

and the metric function
$\psi_{ij}\equiv\psi(d_{ij})=e^{-d_{ij}^{2}/2\sigma^{2}}$ is chosen as a
Gaussian function with a tunable parameter $\sigma$ as it should be a falling
function with increasing the distance between events. $T$ compares the average
distances of pairs of events belonging to two samples of opposite flavour. The
normalisation factor removes the impact of global asymmetries. The distance
between two points in phase space is given by
$d_{ij}=(m_{12}^{2,j}-m_{12}^{2,i},m_{23}^{2,j}-m_{23}^{2,i},m_{13}^{2,j}-m_{13}^{2,i})$,
where the $1,2,3$ subscripts indicate the final-state particles. For no-CPV,
$T$ is expected to be zero, and larger than zero in case of the CPV. This
unbinned technique calculates a $p$-value under the hypothesis of $CP$
symmetry by comparing the nominal $T$ value observed in data to a distribution
of $T$ values obtained from permutation samples, where the flavour of the
${{\cal D}^{0}}$ is randomly reassigned to simulate samples without $CP$
violation. The $p$-value for the no CPV hypothesis is obtained as the fraction
of permutation $T$ values greater than the nominal $T$ value. The $p$-value
from the fitted $T$ distribution can be calculated as the fraction of the
integral of the function above the nominal $T$ value. This approach is used
for the nominal result. If large $CP$ violation is observed, the observed $T$
value is likely to lie outside the range of permutation $T$ values. In this
case the permutation $T$ distribution can be fitted with a generalised extreme
value function

$\displaystyle
f(T;\mu,\delta,\xi)=N\left[1+\xi\left(\frac{T-\mu}{\delta}\right)\right]^{(-1/\xi)-1}$
$\displaystyle\times\exp\left\\{-\left[1+\xi\left(\frac{T-\mu}{\delta}\right)\right]^{-1/\xi}\right\\},$
(15)

with normalisation $N$, location parameter $\mu$, scale parameter $\delta$,
and shape parameter $\xi$. The $p$-value from the fitted $T$ distribution can
be calculated as the fraction of the integral of the function above the
nominal $T$ value. This approach is used for the sensitivity studies. Using
100 permutations, the sensitivity studies are reported in Table 1.

Table 1: Overview of sensitivities to various $CP$ violation scenarios. $\Delta A$ and $\Delta\phi$ denote, respectively, change in amplitude and phase of the resonance $R$. $R$ ($\Delta A$, $\Delta\phi$) |  $p$-value (fit) | Upper limit
---|---|---
$\rho^{0}$ $(4\%$, $0^{\circ})$ | $3.3^{+1.1}_{-3.3}\times 10^{-4}$ | $4.6\times 10^{-4}$
$\rho^{0}$ $(0\%$, $3^{\circ})$ | $1.5^{+1.7}_{-1.4}\times 10^{-3}$ | $3.8\times 10^{-3}$
$\rho^{+}$ $(2\%$, $0^{\circ})$ | $5.0^{+8.8}_{-3.8}\times 10^{-6}$ | $1.8\times 10^{-5}$
$\rho^{+}$ $(0\%$, $1^{\circ})$ | $6.3^{+5.5}_{-3.3}\times 10^{-4}$ | $1.4\times 10^{-3}$
$\rho^{-}$ $(2\%$, $0^{\circ})$ | $2.0^{+1.3}_{-0.9}\times 10^{-3}$ | $3.9\times 10^{-3}$
$\rho^{-}$ $(0\%$, $1.5^{\circ})$ | $8.9^{+22}_{-6.7}\times 10^{-7}$ | $4.2\times 10^{-6}$

By counting the fraction of permutations with a $T$ value above the nominal
$T$ value in the data, a $p$-value of $(2.6\pm 0.5)\times 10^{-2}$ is
extracted. This result is based on 1000 permutations. The results correspond
to a data sample of about 2/fb collected during 2012. Varying the metric
parameter results in the $p$-values listed in Table 2; all the $p$-values are
at the $10^{-2}$ level.

Table 2: Results for various metric parameter values. The $p$-values are obtained with the counting method. $\sigma$ [$\mathrm{GeV}^{2}/\mathrm{c}^{4}$] |  $p$-value
---|---
0.2 | $(4.6\pm 0.6)\times 10^{-2}$
0.3 | $(2.6\pm 0.5)\times 10^{-2}$
0.4 | $(1.7\pm 0.4)\times 10^{-2}$
0.5 | $(2.1\pm 0.5)\times 10^{-2}$

The data sample has been split according to various criteria to test the
stability of the results. Analyses of sub-samples with opposite magnet
polarity, with different trigger configurations, and with fiducial selection
requirements removing areas of high local asymmetry of the tagging soft pion
from the $D^{*}+$ decay all provide consistent results. Various checks have
been performed to ensure there are no asymmetries arising form background
events or detector related asymmetries.

The analysis has the words best sensitivity from a single experiment to local
CPV in this decay.

##  5 Conclusions

LHCb has performed world leading precision measurements in the charm sector.
The searches for direct CPV in two- and multi body decays are consistent with
CP conservation, in agreement with the SM at the current level of precision.

## References

  * [1] LHCb collaboration, R. Aaij et al., Measurement of $CP$ asymmetry in $D^{0}\rightarrow K^{-}K^{+}$ and $D^{0}\rightarrow\pi^{-}\pi^{+}$ decays, JHEP 07 (2014) 041, arXiv:1405.2797
  * [2] LHCb collaboration, A search for time-integrated $C\\!P$ violation in $D^{0}\rightarrow K^{-}K^{+}$ and $D^{0}\rightarrow\pi^{-}\pi^{+}$ decays, LHCb-CONF-2013-003
  * [3] LHCb collaboration, R. Aaij et al., Search for $CP$ violation in $D^{\pm}\rightarrow K^{0}_{S}K^{\pm}$ and $D_{s}^{\pm}\rightarrow K^{0}_{S}\pi^{\pm}$ decays, JHEP 1410 (2014) 25, arXiv:1406.2624
  * [4] M. Williams, Observing CP Violation in Many-Body Decays, Phys. Rev. D84 (2011) 054015, arXiv:1105.5338
  * [5] BaBar Collaboration, B. Aubert et al., Search for CP Violation in Neutral D Meson Cabibbo-suppressed Three-body Decays, Phys. Rev. D78 (2008) 051102, arXiv:0802.4035
  * [6] LHCb collaboration, R. Aaij et al., Search for CP violation in $D^{0}\rightarrow\pi^{-}\pi^{+}\pi^{0}$ decays with the energy test, arXiv:1410.4170, submitted to PLB