Multivariate L\'evy-type drift change detection and mortality modeling

In this paper we give a solution to the quickest drift change detection problem for a multivariate L\'evy process consisting of both continuous (Gaussian) and jump components in the Bayesian approach. We do it for a general 0-modified continuous prior distribution of the change point. Classically, our criterion of optimality is based on a probability of false alarm and an expected delay of the detection, which is then reformulated in terms of a posterior probability of the change point. We find a generator of the posterior probability, which in case of general prior distribution is inhomogeneous in time. The main solving technique uses the optimal stopping theory and is based on solving a certain free-boundary problem. We also construct a Generelized Shiryaev-Roberts statistic, which can be used for applications. The paper is supplemented by two examples, one of which is further used to analyze Polish life tables (after proper calibration) and detect the drift change in the correlated force of mortality of men and women jointly.


INTRODUCTION
Quickest detection problems, often called also disorder problems, arise in various fields of applications of mathematics, such as finance, engineering or economics.All of them address a question how to detect some changes in observed system in an optimal way using statistical methods.One of the main methods was based on the drift change detection using Bayesian approach; see e.g.Shiryaev [25,26], where Brownian motion with linear drift was considered and the drift has been changing according to an exponential distribution.The original problem was reformulated in terms of a free-boundary problem and solved using optimal stopping methods.All details of this analysis are also given in surveys [31,33] (see also references therein).Apart from Baysian method, the minimax approach have also been used.This method is based on identifying the optimal detection time based on so-called cumulative sums (CUSUM) strategy; see e.g.Page [16], Beibel [3], Shiryaev [28] or Moustakides [14] in the Wiener case, or El Karoui et al. [10] in the Poisson case.Many of these quickest detection problems and used methods are gathered in the book of Poor and Hadjiliadis [22].In this paper we choose the first approach.Our first main goal is to perform the analysis of the quickest drift change detection problem for multivariate processes, taking into account the dependence between components.We also allow a general 0-modified continuous prior distribution of the change point.Most of works on the detection problems in Bayesian setting has been devoted to the one-dimensional processes consisting of only continuous (Gaussian) part or only jumps; see e.g.Beibel [2], Shiryaev [26] or [32,Chap. 4] or Poor and Hadjiliadis [22].Only some particular cases of jump models without diffusion component have been already analysed, e.g. by Gal'chuk and Rozovskii [7], Peskir and Shiryaev [19] or Bayraktar et al. [4] for the Poisson process, by Gapeev [8] for the compound Poisson process with the exponential jumps or by Dayanik and Sezer [5] for more general compound Poisson problem.Later, Krawiec et al. [11] allowed observed process to have, apart from diffusion ingredient, jumps as well.This is very important in many applications appearing in actuarial science, finance, etc.Still, all of these results concern one-dimensional case only.This paper removes this limitation.In addition, we assume that a drift change point θ has a general 0-modified continuous prior distribution G.In most works it has been assumed that θ can have only (0-modified) exponential distribution.Such assumption makes the free-boundary problem time-homogenous due to lack of memory property, which is not true in the general case.The methodology used in this paper is based on transferring the detection problem to a certain freeboundary problem.More formally, in this paper we consider the process X = (X t ) t≥0 with (1) where X ∞ = (X ∞ t ) t≥0 and X 0 = (X 0 t ) t≥0 are both independent jump-diffusion processes taking values in R d .We assume that X ∞ and X 0 are related with each other via the exponential change of measure described e.g. in Palmowski and Rolski [17].This change of measure can be seen as a form of the drift change between X ∞ and X 0 with additional change in jump distribution.Later we will see the parameter r, which corresponds to the rate (direction) of disorder that can be observed after time θ.Let θ has an atom at zero with mass x > 0. We choose the classical optimality criterion based on both probability of a false alarm and a mean delay time.That is, in this paper, we are going to find an optimal detection rule τ * ∈ T for which the following infimum is attained where T is the family of stopping times and c > 0 is fixed number.Measure P G will be formally introduced later.Firstly, we transfer above detection problem into the following optimal stopping problem for the a posteriori probability process π = (π t ) t≥0 that also will be formally introduced later.The subscript x associated with E G indicates the starting position of process π equal to x.In the next step, using the change of measure technique and stochastic calculus, we can identify the infinitesimal generator of the Markov process π.This part contains results of independent interest on properties of the posterior process π, that are related to the multidimensionality of the process X.In the classical case with exponential distribution G, π is time-homogenous with generator A. Finally, we formulate the free-boundary value problem, which in the time-homogenous case is as follows with the boundary conditions for some optimal level A * which allows to identify the threshold optimal alarm rule as We first generalize above free-boundary problem and then solve it for two basic models: two-dimensional Brownian motion and two-dimensional Brownian motion with downward exponential jumps.
Our second main goal is to apply the solution of above multivariate detection problem to the analysis of correlated change of drift in force of mortality of men and women.The life expectancies for men and women are widely recognized as dependent on each other.For example, married people live statistically longer than single ones.Since many insurance products are engineered for marriages or couples it is crucial to detect the change of mortality rate of marriages.Indeed, the observed improvements of longevity produce challenges related with the capital requirements that has to be constituted to face this long-term risk and with creating new ways to cross-hedge or to transfer part of the longevity risk to reinsurers or to financial markets.To do this we need to perform accurate longevity projections and hence to predict the change of the drift observed in prospective life tables (national or the specific ones used in insurance companies).In this paper we analyze the Polish life tables for both men and women jointly.We proceed as follows.We take logarithm of the force of mortality of men and women creating a two-dimensional process, modeled then by a jump-diffusion process.This process consists of observed two-dimensional drift that can be calibrated from the historical data and a random zero-mean Lévy-type perturbation.Based on previous theoretical work we construct a statistical and numerical procedure based on the generalized version of the Shiryaev-Roberts statistic introduced by Shiryaev [25,26] and Roberts [23], see also Polunchenko and Tartakovsky [20], Shiryaev [29], Pollak and Tartakovsky [21] and Moustakides et al. [15].Precisely, we start from a continuous statistic derived from the solution of the optimal detection problem in continuous time.Then we take discrete moments 0 < t 1 < t 2 < . . .< t N , construct an auxiliary statistic and raise the alarm when it exceeds certain threshold A * identified in the first part of the paper.The set-up used in examples is, however, simplified compared to the theory presented in the previous sections.Applications focus mainly on multi-dimensionality of presented problem, to see how one can analyse mortality of men and women jointly.The distribution of change time θ is limited to classical (0-modified) exponential.The paper is organized as follows.In Section 2 we describe basic setting of the problem, introduce main definitions and notation.In this section we also formulate main theoretical results of the paper.Section 3 is devoted to the construction of the Generalized Shiryaev-Roberts statistic.To apply it, we first need to find some density processes related to the processes X prior and post the drift change.This is done in Section 3 as well.Particular examples are analyzed in Section 4. Next, in Section 5, we give an application of the theoretical results to a real data from life tables.We finish our paper with some technical proofs given in Section 6.

MODEL DESCRIPTION AND MAIN RESULTS
The main observable process is a regime-switching d-dimensional process X = (X t ) t≥0 .It changes its behavior at a random moment θ in the following way: (2) where X ∞ and X 0 are two different independent Lévy processes related with each other via exponential change of measure specified later.The random time θ is independent of the pre-and post-change random processes.
We assume the following model when the post-change drift equals r.The process that we observe after the change of drift is a d-dimensional processes X 0 = (X 0 t ) t≥0 = (X 0 t,1 , . . ., X 0 t,d ) t≥0 defined as (3) where T is a vector of standard independent Brownian motions, • σ = (σ i,j ) i,j=1,...,d is a matrix of real numbers, responsible for the correlation of the diffusion components of X 0 t,1 , . . ., X 0 t,d , we assume that σ ii > 0 for all i = 1, . . ., d, • r = (r 1 , . . ., r d ) T is a vector of an additional drift, • N 0 = (N 0 t ) t≥0 is a Poisson process with intensity µ 0 , • (J 0 k ) k≥1 is a sequence of i.i.d.random vectors responsible for jump sizes; we denote each coordinate of J 0 k by J 0 k,i for i = 1, . . ., d and its distribution by F 0 i with mean m 0 i ; we also denote by F 0 a joint distribution of vector J 0 k and by m 0 = (m 0 1 , . . ., m 0 d ) T its mean.
We assume that all components of X 0 t are stochastically independent, i.e.W 0 t , N 0 t and the sequence (J 0 k ) k=1,2,... are independent.
Similarly, we assume that the process that we observe prior the drift change is a d-dimensional process where t≥0 is a vector of standard independent Brownian motions, • matrix σ is the same as for the process X 0 , We denote W t = W ∞ t∧θ + W 0 (t−θ) + which is a Brownian motion as well.To formally construct the model with a drift change described above, we follow the ideas of Zhitlukhin and Shiryaev [37].Precisely, we consider a filtered measurable space (Ω, F, {F t } t≥0 ) with a rightcontinuous filtration {F t } t≥0 , on which we define a stochastic system with disorder as follows.First, on a probability space (Ω, F, {F t } t≥0 ) we introduce two probability measures P ∞ and P 0 with their restrictions to F t given by P ∞ t := P ∞ | Ft and P 0 t := P 0 | Ft .We assume that for each t ≥ 0 the restrictions P ∞ t and P 0 t are equivalent.The measure P ∞ corresponds to the case when there is no drift change in the system at all and P 0 describes the measure under which there is a drift r present from the beginning (i.e. from t = 0).In the following we assume that both measures correspond to laws of the processes X ∞ and X 0 described above, respectively.We also introduce a probability measure P that dominates P ∞ and P 0 and such that the restriction P t := P| Ft is equivalent to P ∞ t and P 0 t for each t ≥ 0. We define the Radon-Nikodym derivatives (5) Finally, for any fixed s ∈ (0, ∞), taking a consistent family of probability measures (P by the Kolmogorov's existence theorem we can define measures P (s) such that t .Note that for t < s the following equality holds since disorder after time t does not affect the behavior of the system before time t.We consider Bayesian framework, that is, we assume that the moment of disorder is a random variable θ with a given distribution function denoted by G(s) on (R + , B(R + )).We assume that G(s) is continuous for s > 0 with right derivative G (0) > 0. We define all quantities on an extended filtered probability space (Ω, F, {F t } t≥0 , P G ) such that ( 7) Observe that measure P G describes formally the process X defined in (2).
In the problem of the quickest detection we are looking for an optimal stopping time τ * that minimizes certain optimality criterion.We consider a classical criterion, which incorporates both the probability of false alarm and the mean delay time.Let T denote the class of all stopping times with respect to the filtration {F t } t≥0 .Our problem can be stated as follows: Problem 1.For a given c > 0 calculate the optimal value function and find the optimal stopping time τ * for which above infimum is attained.
Above E G means the expectation with respect to P G .The key role in solving this problem plays a posterior probability process π = (π t ) t≥0 defined as ( 9) We denote x := π 0 = G(0) and add a subscript x to E G x to emphasize it.Using this posterior probability, one can reformulate criterion (8) into the following, equivalent form: Problem 2. For a given c > 0 find the optimal value function and the optimal stopping time τ * such that That is, formally, the following result holds true.

Lemma 1. The criterion given in Problem 1 is equivalent to the criterion given in Problem 2.
Although the proof follows classical arguments, we added it in Section 6 for completeness.Below we formulate the main theorem that connects Problem 2 to the particular free-boundary problem.It is based on the general optimal stopping theory in the similar way as Theorem 1 in Krawiec et al. [11], which it extends.However, for the general (continuous for s > 0 with right derivative G (0) > 0) distribution G(s) of the moment θ, the optimal stopping problem and its solution are time-dependent.The problem reduces to time-independent case for the (0-modified) exponential distribution G.We will prove it in Section 6.

Theorem 1. Let ∂
∂t + A be a Dynkin generator of the Markov process (t, π t ) t≥0 .Then the optimal value function V * (x) from the Problem 2 equals f 0 (x), where f t (x) solves the free-boundary problem with the boundary conditions Furthermore, the optimal stopping time for the Problem 2 is given by If G is the (0-modified) exponential distribution, then V * (x) solves above free-boundary problem for the unique point A * ∈ (0, 1] not depending on time with the optimal stopping time given by Further, in this case f t (x) = f 0 (x) = f (x) and additionally the following condition holds See also Peskir  It is known that the Dynkin generator is an extension of an infinitesimal generator in the sense of their domains.Following [18, Chap.III] and discussion done on page 131 of [18] (see also the proof of [9, Prop.2.6]) we can conclude that the optimal value function V * (t, x) satisfies ( 10) where ∂ ∂t + A is an infinitesimal generator as long as there exists unique solution of (10) lying in the domain of infinitesimal generator.Now to formulate properly above free-boundary problem, we have to identify the infinitesimal generator ∂ ∂t + A and its domain.They are given in next theorem.We use notation Theorem 2. The infinitesimal generator of the Markov process (t, π t ) t≥0 is given by ∂ ∂t f t (x) + Af t (x) for and for functions We will prove this theorem later in Section 6. Assume that we can find unique solution of ( 10)- (12) in the class C 2 with A given in (16) then by above considerations it follows that this solution equals the value function V * (t, x).Therefore in the final step we focus on the simple time-homogeneous case of exponential time change case, then we solve uniquely (10)- (12) for some specific choice of model parameters, and finally, find the optimal threshold A * and hence the optimal alarm time.This allows us to construct a Generalized Shiryaev-Roberts statistic in this general set-up.Later we apply it to detect the changes of drift in joint (correlated) mortality of men and women based on life tables.
We will give another representation of the process π in terms of the process L = (L t ) t≥0 defined by ( 18) To find above Radon-Nikodym derivative such that process X defined in (2) indeed admits representation (3) under the measure P 0 and (4) under P ∞ , we assume that for given r = (r 1 , . . ., r d ) ∈ R d the following relation holds where for x = (x 1 , . . ., x d ) the function h r (x) : R d → R is given by Above the coefficients z r,1 . . .z r,d solve the following system of equations: Theorem 3. Assume that (19) holds for given r ∈ R d and that the Radon-Nikodym derivative L = (L t ) t≥0 defined in ( 18) is given by Then the process X defined in (2) admits representation (3) under the measure P 0 and (4) under P ∞ .
The proof will be given in Section 6.
Having the density process L defined in (18) identified in above theorem, we introduce an auxiliary process (24) Then by ( 6), ( 17) and ( 18) the following representation of π t holds true (25) , where the last equality follows from the definition of L (s) t in (6) for t < s.By the It ô's formula applied to (24) we obtain that ψ 0 t solves the following SDE ( 26) The construction of the classical Shiryaev-Roberts statistic (SR) is in detail described and analyzed e.g. by Shiryaev [29], Pollak and Tartakovsky [21] and Moustakides et al. [15].In this paper we consider Generalized Shiryaev-Roberts statistic (GSR).We start the whole construction from taking the discrete-time data X ti ∈ R d observed in moments 0 = t 0 < t 1 < . . .< t n , where n is a fixed integer.We assume that t i − t Considering a discrete analogue of ( 24) we define the following statistic where from equation ( 22) we take for n > 0 and L 0 = 1.Above G(0) = x corresponds to an atom at 0. For convenience it can be also calculated recursively as follows: Recall from Theorem 1 that the optimal stopping time is given by for some optimal level A * .Therefore from identity (25) we can introduce the following Generalized Shiryaev-Roberts statistic and raise the alarm of the drift change at the optimal time of the form We emphasize that the GSR statistic is more appropriate in longevity modeling analyzed in this article than the standard one (i.e.SR).Indeed, the classical statistic is a particular case when θ has an exponential distribution with parameter λ tending to 0. The latter case corresponds to passing with mean value of the change point θ to ∞ and hence it becomes conditionally uniform, see e.g.Shiryaev [29].Still, in longevity modeling it is more likely that life tables will need to be revised more often and therefore keeping dependence on λ > 0 in our statistic seems to be much more appropriate.For the similar reasons we also prefer to fix average moment of drift change θ instead of fixing the expected moment of the revision time τ .To apply above strategy we will focus on the exponential time of drift change.In this case we have to identify the optimal alarm level A * in the first step and hence we have to solve uniquely the freeboundary value problem ( 10) - (12).We analyze two particular examples in the next section.

EXAMPLES
4.1.Two-dimensional Brownian motion.Consider the process X without jumps (i.e. with jump intensities µ ∞ = µ 0 = 0).In terms of processes X 0 and X ∞ given in (3) and ( 4) it means that (27) Assume that Then the first coordinate X 0 t,1 is a Brownian motion with drift and with variance σ 2 1 and the second coordinate X 0 t,2 is also a Brownian motion with drift and with variance σ 2 2 .The correlation of the Brownian motions on both coordinates is equal to ρ. Process X ∞ has similar characteristics but without any drift.Next, assume that, conditioned on θ > 0, θ is exponentially distributed with parameter λ > 0, i.e.
Then the generator of process π according to ( 16) is equal to where z r,1 and z r,2 solve the following system Our goal is to solve the boundary value problem ( 10) -( 12) where generator A is given by (28).Note that the system (10) takes now the following form where Observe that above equations allow us to refer to the classical Shiryaev problem, with our constant B included.Hence, from Shiryaev [25,31] it follows that solution of above equation is given by where The exact values of function y(x) can be found numerically, while the threshold A * can be found from the equation y(A * ) = −1, which is the boundary condition (12).

Two-dimensional Brownian motion with one-sided jumps.
The second example concerns similar 2-dimensional Brownian motion model as in the previous example, but with additional exponential jumps.Assume that µ ∞ , µ 0 > 0 and (29) w j e −yj /wj I(y j ≥ 0)dy.
In other words, jump sizes on each coordinate j ∈ {1, 2} of process X ∞ are independent of each other and distributed exponentially with mean w j > 0. Additionally, we assume as in the previous example that Jump distribution given by ( 29) together with Theorem 3 allows us to formulate the following lemma.
Lemma 2. Assume that jump distribution F ∞ of the process X ∞ is given by ( 29) and jump intensity is equal to µ ∞ .Assume also that there exists a vector z r = (z r,1 , . . ., z r,d ) satisfying the system (21) such that (∀ 1≤j≤2 )(|w j z r,j | < 1).Then the following distribution function F r and intensity µ 0 satisfy the condition (19): 1 − w j z r,j w j e −yj /( w j 1−w j z r,j ) I(y j ≥ 0)dy, Proof.By the combination of ( 19) and ( 20) we obtain that −yj /( w j 1−w j z r,j ) I(y j ≥ 0)dy, which can be rearranged to 1 − w j z r,j w j e −yj /( w j 1−w j z r,j ) I(y j ≥ 0)dy.Now it is sufficient to observe that above formula is equal to the product µ 0 F r given by (30) and that F r is indeed a proper distribution by the assumption that (∀ 1≤j≤2 ) (|w j z r,j | < 1).

Remark 1.
Considering the jump distributions F ∞ and F r given by ( 29) and (30), the system (21) consists of equations .
Remark 2. Distribution F r given by (30) has similar characteristics to F ∞ .More precisely: jumps on both coordinates X 0 t,1 and X 0 t,2 are independent, exponentially distributed with means The generator A given by ( 16) for jump distributions specified above can be expressed as The integral part of A can be further simplified.For α 1 , α 2 > 0 we define the following integrals α j e −αj yj dy and α j e αj yj dy.

Remark 3.
In Lemma 3 we restrict calculations to the case β 1 = β 2 , but similar transformations of the integral may be made for the case β 1 = β 2 as well.The difference will appear in the distribution of random variable S 0 present in the proof (see Section 6) being the sum of two exponential random variables.
Denote γ i := 1 zr,iwi for i ∈ {1, 2}.Then from Lemma 3 the generator A given in ( 31) can be rewritten as follows Equation Af (x) = −cx in the free-boundary value problem can be further simplified to get rid of the integrals and then solved numerically to find the threshold A * .We believe that this particular case may be finally solved numerically in a similar way as in the numerical analysis described in Krawiec et al. [11], since here we obtain equation of the same order and similar characteristics.However, in this article we focus our applications on the previous example, which is used in practice in the next section.
Remark 4. The results of above example are derived under the assumption of positive exponential jumps.However, the whole analysis can be also conducted for negative exponential jumps, i.e. for the distribution e yj /wj I(y j ≤ 0)dy.
Then we can use part of Lemma 3 concerning I 0 − (x) to derive the generator A given by

APPLICATION TO THE FORCE OF MORTALITY
Now we are going to give an important example of applications, which concerns modeling of the force of mortality process.We will analyze the joint force of mortality for both men and women.We observe this process over the past decades and check if and when there have been significant changes of drift.To achieve this goal, we introduce two-dimensional process of the force of mortality µ := (µ t ) t≥0 = ((µ 1 t , µ 2 t )) t≥0 .We interpret this process as follows: • the first coordinate µ 1 t represents force of mortality of men, while the second one µ 2 t represents force of mortality of women (of course they are correlated), • the time t runs through consecutive years of life tables, e.g. if t = 0 corresponds to the year 1990, then t = 10 corresponds to the year 2000, • the age of people is fixed for a given process µ, i.e. if µ 0 concerns 50-year old men and women, then µ 10 also concerns 50-year old men and women, but in another year.
The representation of the force of mortality process is given by ( 32) where log μt := (log µ 1 t , log µ 2 t ) is a deterministic part equal to log μt = a 0 + a 1 t.
Above a 0 = (a 1 0 , a 2 0 ) is a known initial force of mortality vector of men and women and a 1 = (a 1 1 , a 2 1 ) is a vector of a historical drift per one year.It is worth to mention here that our model is similar to the Lee-Carter model (for fixed age ω, cf.[13] ): where a ω is a chosen number, k t is certain univariate time series and ω,t is a random error.However, Lee-Carter method focuses on modelling the deterministic part of the force of mortality, while our detection procedure concerns controlling the random perturbation in time, precisely the moment when it substantially changes.This model is univariate as well in contrast to our two-dimensional mortality process.In our numerical analysis the stochastic part X t will be modeled by the two-dimensional Brownian motion analyzed in Example 4.1.We apply this model to the life tables downloaded from the Statistics Poland website [35].The first step concerns the model calibration.We start with some historical values of the force of mortality μ0 , . . ., μn , where each μi = (μ i,1 , μi,2 ) is a two-dimensional vector (one coordinate for women and one for men).We estimate a 1 as a mean value of log-increments of μ0 , . . ., μn .Precisely, where A little more attention is needed to calibrate the stochastic part X, which includes correlation.Denote (33) Xi = log μi − a 0 − â1 i, i = 0, . . ., n and the increments (34) x i := Xi+1 − Xi , i = 1, . . ., n.
There are still some model parameters that have to be chosen a priori.In particular, we have to declare the anticipated incoming drift r, the probability x = π 0 = P G (θ = 0) that the drift change occurs immediately, the parameter λ > 0 of the exponential distribution of θ and parameter c present in criterion stated in the Problem 2.
We assume their values at the following level: • λ = 0.1.It is the reciprocal of the mean value of θ distribution conditioned to be strictly positive.Such choice reflects the expectation that the drift will change in 10 years on average.
• x(= P G (θ = 0)) = 0.1.This parameter should be rather small (unless we expect the change of drift very quickly).• c = 0.1.It is the weight of the mean delay time inside the optimality criterion stated in the Problem 1.It reflects how large delay we can accept comparing to the risk of false alarm.We have chosen rather small value and connected it to λ by choosing c = λ.• Drift incoming after the moment θ -we have connected the anticipated value of r to σ by r = (σ 1 , σ 2 ).In practice we suggest to adjust the choice of r to the analysis of sensitivity of e.g.price of an insurance contract.
In Table 1 we sum up all parameters that were used (both calibrated and arbitrary chosen ones) in the numerical analysis.The calibration interval was set to years 1990 -2000.In the Figure 1 we present exemplary plot of the force of mortality for women at age 60 through years 1990 -2017.Most of the time it is decreasing, but we can observe a stabilization period around years 2002 -2009.According to (32) we first take logarithm of the force of mortality, separate deterministic linear part and then model the remaining part by the process X given by (27). Figure 2 presents historical observations of this remaining part for the same data as in the Figure 1.The results of the detection algorithm for the force of mortality of 60-year old men and women jointly are presented in the Figure 3.The change of drift for given parameters was detected in year 2006 (red vertical line in the first two plots).The threshold A * for the optimal stopping time is here equal to 0.85, which is indicated by the red horizontal line in the third plot presenting values of π = (π t ) t∈{1990,...,2017} .Note that calibration of parameters (including historical drift) has been done for interval 1990 -2000, when the force of mortality was mostly decreasing.After year 2002 it stayed at a stable level for several years, which was detected as a change of drift.This change of behavior is even more evident  2, where we can observe that process X is mostly increasing through the years 2002 -2009.This example shows that our detection method does not necessarily rise the alarm after the first observed deviation, but rather after it becomes more evident, that the change of drift actually has happened.Therefore, it copes well with cases of gradually changing drift, as long as eventually observed process significantly deviates from the model.An important note need to be given at the end.This procedure is strongly dependent on parameters chosen to the model -e.g.post-change drift vector r, which was chosen depending on σ 1 and σ 2 , to give appropriate order of magnitude.A full analysis of the impact of individual parameters on the results would significantly extend this article.However, it may be the subject of further research in further articles, developing the applications of the detection method described here.

PROOFS
Proof of Lemma 1.Note that Moreover, observe that by Tonelli's theorem we have: Putting together (35) and ( 36) completes the proof.
Proof of Theorem 1.We start from the observation that process ((s, π s )) s≥0 is Markov which follows from Theorem 2. Let Then V * (x) = V * (0, x).Moreover, for fixed t ≥ 0 the optimal value function V * (t, x) is concave, which follows from [11,Lem. 3] and the assumption that distribution function G(t) of θ is continuous for t > 0. Observe that from Theorem 2 it follows that ((s, π s )) s≥0 is stochastically continuous and and by [18, Chap.III] the optimal value function V * (t, x) satisfies the following system (38) where A is a Dynkin generator.Using the same arguments like in the proofs of [11,Lem. 6 and Lem. 7] we can prove the boundary conditions f Finally, if G is exponential, then (π s ) s≥0 is Markov by Theorem 2. Now we can do the same arguments but taking simply V * (x) instead V * (t, x) above.Moreover, (15) follows from the same arguments like in the proof of [11,Lem. 6].This completes the proof.Proof of Theorem 2. First, we will find the SDE which is satisfied by process π t .By the definition of the process X for each i = 1, . . ., d, we get Denote the continuous part of the process by an additional upper index c.Then For the process L t given in (22), by the It ô's formula we obtain where By (26) we conclude that Recall that by (25) we have .
Then, using It ô's formula once again we obtain Moreover, Together with the system of equations ( 21) it produces Using the It ô's formula one more time completes the proof.
Proof of Theorem 3. The proof is based on the technique of exponential change of measure described in Palmowski and Rolski [17].Firstly, we will prove that the process (L t ) t≥0 satisfies the following representation for the function h(x) := h r (x) given in (20), where A ∞ is an extended generator of the process X under P ∞ and h is in its domain since it is twice continuously differentiable.Then from Theorem 4.2 by Palmowski and Rolski [17] it follows that the generator of X under P 0 is related with A ∞ by (40) On the other hand, from the definition of the infinitesimal generator or using the Theorem 31.5 in Sato [24] it follows that for twice continuously differentiable function f (x 1 , . . ., x d ) : R d → R generators A ∞ and A 0 are given by For h r (x) given by ( 20) we obtain h r (X t ) h r (X 0 ) = exp z r,j (X t,j − X 0,j ) Further, since ∂h r ∂x i = z r,i h r and Hence, we obtain z r,j (X t,j − X 0,j ) − K r t    and thus L 0 t given in (22) indeed satisfies the representation (39) for function h r (x) given by (20).
To finish the proof it is sufficient to show that the generator A 0 given by (42) indeed coincides with the generator given in (40) for h(x) = h r (x).First, by (23) we get Second, (41) produces ∂ 2 f ∂x i ∂x j + ∂f ∂x i z r,j + ∂f ∂x j z r,i + f z r,j z r,i (σσ T ) i,j Hence (f (x + y) − f (x))µ 0 F 0 (dy).
FIGURE 1. Force of mortality of women aged 60 in 1990-2017