### Växjö

### University

School of Mathematics and System Engineering
**Reports from MSI - Rapporter från MSI**

### Hatice Yalman

### Change Point Estimation for Stochastic

### Differential Equations

Bachelor’s Thesis Mathematics

2009

**Abstract**

A stochastic differential equation driven by a Brownian motion where the dispersion is determined by a parameter is considered. The parameter undergoes a change at a certain time point. Estimates of the time change point and the parameter, before and after that time, is considered.The estimates were presented in Lacus 2008. Two cases are consid-ered: (1) the drift is known, (2) the drift is unknown and the dispersion space-independent. Applications to Dow-Jones index 1971-1974 and Goldmann-Sachs closings 2005– May 2009 are given.

**Key-words:** Brownian motion, stochastic differential equations, Ornstein-Uhlenbeck,
change points, estimates, simulations, closings, returns, Dow-Jones, Goldman-Sachs.

**Abstract**

Bu tezi yazarken amacimiz stochastic differensiyel esitliginden yararlanarak parametreye bagli olan gercek degisimin oldugu degeri bulmaktir. Bu degeri bulmak icin once bazi verilerin bilindigi, sonra bu degerlerin bilinmedigi kabul edilerek gercek degisimin oldugu degeri bulmak icin bazi islemler yapilmistir.

**Key-words:**

**Contents**

**1** **Overview** **1**

**2** **Introduction to the Brownian motion and stochastic differential equations** **1**

**3** **Change point estimation** **1**

3.1 Change point estimation with known drift . . . 2 3.2 Change Point Estimation With Unknown Drift . . . 8

**4** **appendix** **19**

4.1 Typical simulation of an SDE with change-point . . . 19 4.2 Known drift . . . 19 4.3 Unknown drift . . . 20

**1**

**Overview**

Much of this work is inspired by [2]

In section 2;We gave a brief introduction of a Brownian motion and Stochastic differ-ential equation.

In section 3;The model and the estimates first with known drift and secondly with unknown drift.The estimates are applied to real data.

**2**

**Introduction to the Brownian motion and stochastic differential**

**equations**

*A random variable is said to be N(0,1) if it has density function f(x) = e−x*2/2/√2π,

*x*_{∈ R.}

*The process W* _{= {W}t*: t* >= 0} is a Brownian motion if

*(i) Wt* *is continuous and W*0= 0

*(ii) the value of Wt*/√*t is N(0,1)*

*(iii) the increment Ws+t−Ws* is N*(0,t), i.e. (Ws+t−Ws*)/√*t is N(0,1), and the *

incre-ments are independent of F*s, the history of what the process did up to time s.*

*These are both necessary and sufficient conditions for a process W to be a Brownian*
motion.

Shortly: a Brownian motion is a continuous-time stochastic process_{{W}t*: t*≥ 0} such

*that Wt*/√*t is N(O,1).*

A SDE is a differential equation typically written in a form

*dXt= b(x)dt +*σ*(x)dW* (2.1)

*where W is a BM. We can do an infintesimal interpretation: As the time goes from t to*

*dt,X is moving a time dX given by (2.1), where*

*dW* =√*dtN*(0, 1). (2.2)

*or in other words, W _{t}_{+dt}*

*= Wt+ dW , where dW is given by (2.2) and similarly, Xt+dt*=

*Xt+ dX, where dX is given from (2.1). There are number of different definitions of what*

a solution to (2.1) using stochastic integrals but we don‘t need that for this bachelor thesis.

**3**

**Change point estimation**

In this thesis we study estimation of a volatility parameter θ for stochastic differential equations which changes value after a distinct time point.

The crucial thing is that the parameterθ shifts value. In factθ =θ*(t) depends on time.*

In this thesis, for estethical reasons we we consequently writeθ instead ofθ*(t). Before a*

timeτ0it has a valueθ1and afterτ0,θ =θ2:

θ =

θ1*, t ≤*τ0

θ2*, t >*τ0

(3.1) We assume throughout, without loss of generality, thatθ1 andθ2 are non-negative. The

reason is thatσ*(Xt)dW has the same distribution as −*σ*(Xt)dW (talking infinitesmal; the*

0 0.2 0.4 0.6 0.8 1 −4 −3 −2 −1 0 1 2 3 4 5

*Figure 3.1: Simulated data X given by dXt= Xtdt*+θ*dWt*, whereθ *= 1 for t <*τ0= 0.5

andθ * _{= 10 for t ≥}*τ0

**3.1** **Change point estimation with known drift**

In this section we consider estimation of the parameterθ and τ0 given by (3.1) for the

SDE

*dXt= b(Xt)dt +*θσ*(Xt)dWt*, (3.2)

*where the drift b and the dispersion*σ are known functions.

In this chapter we produce approximate solutions of (3.2) with Euler-Maruyama
sim-ulations during a time interval*[0, T ]*

*Xi*+1*= Xi+ b(Xi*)∆*n*+θσ*(Xi)(Wi*+1*−Wi*) (3.3)

where∆*i= ti*+1*− ti= T /(n − 1), Wi*+1*−Wi*is a normal (Gaussian) random number with

mean zero and standard deviation√∆*i. By the notation Xiand Wiin (3.3 we mean Xti*and

*Wti*, respectively. That notation was also made due to estethical reasons. A simulation of

(3.3) is in Figure 3.1. The reason why Euler-Maruyama simulations are performed is to check the validity of estimates ofθ1,θ2andτ0presented below.

*The problem setup is now as follows. Assume we know the functions b(x) and*σ*(x),*

*initial value X*0and a trajectory*{Xt} at gridpoint 0 = t*0*< t*1*< . . . < tn= T , but we don’t*

knowθ1,θ2, and change time pointτ0. How do we then estimateθ1,θ2andτ0?

To estimateθ1, θ2 andτ0 we do as follows. First we explain the algorithm and then

summarize it. Let

*Zi*=
*(Xi*+1*− Xi) − b(Xi*)∆*n*
√_{∆}
*n*σ*(Xi*)
(3.4)
Note that
*Zi*=θ
*Wi*+1*−Wi*
√_{∆}
*n*
,

*which, since Wi*+1*−Wi*is √∆*nN(0,1)-distributed, implies that Zi*is θN(0,1)-distributed.

*See Figure 3.2 for Zi*-values from (3.2).

0 0.2 0.4 0.6 0.8 1 −30 −20 −10 0 10 20 30

*Figure 3.2: Zigiven by (3.4) with data Xi*as simulated in Figure 3.1.θ *= 1 for t <*τ0and

θ *= 10 for t ≥*τ0= 0.5
0 0.2 0.4 0.6 0.8 1
0
5
10
15
20
25
30

Figure 3.3:* _{|Z}i| given by (3.4) with data Xi*as simulated in Figure 3.1.θ

*= 1 for t <*τ0and

θ *= 10 for t ≥*τ0*= 0.5. |Zi*|2is in mean equal toθ2so*|Zi*| is aboutθ but it is well-known

0 0.2 0.4 0.6 0.8 1 0 10 20 30 40 50 60

*Figure 3.4: Skgiven by (3.5) with data Xi*as simulated in Figure 3.1.θ *= 1 for t <*τ0= 0.5

andθ * _{= 10 for t ≥}*τ0.

*k*0 would be known, then an estimate ofθ1would be simply

ˆ
θ1=
s
∑*k*0
1 *Zi*2
*k*0

andθ2would be estimated by

ˆ
θ2=
s
∑*n*
*k*0+1*Z*
2
*i*
*n _{− k}*0
.

That was the case whenτ0*= tk*0 was known. Now, the crucial thing is that we don’t know

τ0. How do we do then?

The algorithm with unknownτ0*is as follows. For k= 1, . . . , n, let*

*Sk*=
*k*

### ∑

*i*=1

*Z _{i}*2. (3.5)

*For k= 1, . . . , n, let ˆ*θ1*(k) be an estimate of*θ1,

ˆ
θ1*(k) =*
s
∑*n*
*i*=1*Zi*2
*k* =
r
1
*kSk*,

and let ˆθ2*(k) be an estimate of*θ2,

*What creature is that? U _{k}*2 is a simultaneous measure of the error between ˆθ1

*(k) and*θ1

and between ˆθ2*(k) and*θ2*. If i≤ k*0*, then Z _{i}*2≈θ

_{1}2

*and if i> k*0

*, then Z*2≈θ

_{i}_{2}2.

*Now we should select that k which gives the smallest error, it means the following:*
*The estimate of the change-point index k*0is

*ˆk*0= arg min
*k* {
*k*

### ∑

*i*=1

*(Z*2

_{i}_{−}θ

_{1}2)2+

*n*

### ∑

*i=k+1*

*(Z*2

_{i}_{−}θ

_{2}2)2

_{}}

The estimate of the change time pointτ *is then t _{ˆk}*

0.

*To summarize we have the following algorithm with known drift b and known *
disper-sionσ:
*1. For i= 1, . . . , n, let*
*Zi*=
*(Xi*+1*− Xi) − b(Xi*)∆*n*
√_{∆}
*n*σ*(Xi*)
*2. For k= 1, . . . , n, let*
*S _{k}*=

*k*

### ∑

*i*=1

*Z*2

_{i}*3. For k= 1, . . . , n, let*ˆ θ1

*(k) =*s ∑

*n*

*i*=1

*Zi*2

*k*= r 1

*kSk*, and ˆ θ2

*(k) =*s ∑

*n*

*i=k+1Zi*2

*n*= r 1

_{− k}*n*).

_{− k}(Sn− Sk*4. For k= 1, . . . , n, let*

*U*2=

_{k}*k*

### ∑

*i*=1

*(Z*2

_{i}_{− ˆ}θ1

*(k)*2)2+

*n*

### ∑

*i=k+1*

*(Z*2

_{i}_{− ˆ}θ2

*(k)*2)2, 5.

*ˆk*0= arg min

*k*{

*k*

### ∑

*i*=1

*(Z*2

_{i}_{−}θ

_{1}2)2+

*n*

### ∑

*i=k+1*

*(Z*2

_{i}_{−}θ

_{2}2)2

_{}.}

Figures 3.8 shows the approximate estimates ˆθ1, ˆθ2, ˆτ0. In figures from 3.5-3.9 are graphs

*of U*2, ˆθ1, ˆθ2, ˆτ

There are some asymptotic properties for how good the estimates are for the case
*when b and*σ are known. We will not discuss the facts, just present them. For further
explanations, see [2].

* Fact 3.1 ([2]). Under the null-hypothesis H*0:θ1=θ2, then

*r n*

2*|Dk*|

*d*

*→ |W*0(τ_{)|,}

where* _{{W}*0(τ

_{), 0 ≤}τ

_{≤ 1} is a Brownian bridge ending and starting at zero.}

*By saying ’under the null-hypothesis H*0: θ1=θ2’, we speak in a statistical

0 0.2 0.4 0.6 0.8 1 0.95 1 1.05 1.1 1.15 1.2 1.25x 10 7

*Figure 3.5: U*2 *for X given by dXt* *= Xtdt*+θ*dWt*, where θ *= 1 for t < k*0= 0.5 and

θ *= 10 for t ≥*τ0. Estimated change point ˆτ0*= 0.500. For each time point tk* *the U _{k}*2

-value is computed. 0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7

Figure 3.6: Estimates ofθ1*for X given by dXt= Xtdt*+θ*dWt*, whereθ *= 1 for t <*τ0=

0.5 andθ * _{= 10 for t ≥}*τ0. The estimate is as best at or just before the time change-point

0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12 14

Figure 3.7: Estimates ofθ2*for X given by dXt= Xtdt*+θ*dWt*, whereθ *= 1 for t <*τ0=

0.5 andθ * _{= 10 for t ≥}*τ0. The estimate is as best at or just after the time change-pointτ.

*For each time point tk* there is an estimate ˆθ2*(k).*

0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 12

Figure 3.8: Estimate ofθ *for X given by dXt= Xtdt*+θ*dWt*, whereθ *= 1 for t < k*0= 0.5

and θ * _{= 10 for t ≥}* τ0. Estimated change time point ˆτ0 = 0.500, ˆθ1 = 0.9479, ˆθ2=

10.1200. The estimates ˆθ1and ˆθ2 are drawn with a red (dashed) line. The red (dashed)

line may be difficult to read, due to the precision of the estimates. The red (dashed) line almost coincides with green (solid) line. The exact θ1 and θ2 are drawn with a green

(solid) line. Note that the final estimates ˆθ1 and ˆθ2 are selected at the estimated change

time point ˆτ0. The shaky curve is not needed for the final estimate. It just shows how the

**Fact 3.2 ([2]). The estimator ˆ**τ0*= ˆk*0*/n satisfies*

| ˆτ0−τ0*| = n*−1/2(θ22−θ12)−1*Op*(

p

*log n*).

Moreover; for anyβ _{∈ (0,1),}

*n*β( ˆτ0−τ0)
*P*
→ 0.
Finally
ˆ
τ0−τ*= O(*
1
*n*(θ2
2−θ12)
).
**Fact 3.3 ([2]).**
√
*n*
_{ˆ}
θ2
1−θ12
θ2
1−θ12
*d*
*→ N(0,*Σ)
where
Σ= 2 ˆτ0−1θ04 0
0 2_{(1 −}τ0)−1θ_{0}4
is a covariance matrix.

**3.2** **Change Point Estimation With Unknown Drift**

In this section we consider estimation of the parameterθ for the SDE

*dXt= b(Xt)dt +*θ*dWt*. (3.6)

*Here the drift b is unknown and the dispersion*σ is simply one. In the previous subsection
*the b was assumed to be a known function.*

Note that (3.6) is more general than one may think at the first sight and may be
ap-plied to case (3.2) with non-constantσ*. The equation (3.2) with solution Xt* can, under

non-degeneracity ofθσ (true if for example θσ is bounded away from zero) always be
transformed into one with diffusion coefficient equal to one by applying the Lamperti
transform,
*Yt= F(Xt*) =
Z *X _{t}*

*z*1 θσ

*(u)du*

*Here z is any arbitrary value in the state space of X. Indeed,the process _{{Y }}t* solves the

stochastic differential equation

*dYt* *= bY(Yt)dt + dWt*,
where
*bY(y) =*
*b(F*−1*(y))*
θσ*(F*−1* _{(y))}*−
1
2θσ

*x(F*−1

*which can also write as*

_{(y))}*dYt* = (

*b(Xt*)

θσ*(Xt*)−

1

2θσ*x(Xt))dt + dWt*.
*The process Zt*=θ*Yt* solves the SDE

As to be seen, theθ unfortunately now appears at the drift, which makes it not that easy.
Ifθ2σ*x(Xt*) is small it makes it easier.

One well-known example of (3.6) is the Ornstein-Uhlenbeck-process,

*dX*= (α_{−}β*Xt)dt +*θ*dWt*

whereα,βare constants andθ is a parameter that changes value at some certain timepoint as before. The Ornstein-Uhlenbeck-process is, with other parametrization also called the Vasicek model which is used to describe interest rates.

*b(x) can be estimated nonparametrically with*

*ˆbn(x) =*
∑*n*−1
*i*=1 *K*(
*x _{−X}i*

*hn*

*)(Xi*+1

*− Xi*) ∆

*n*∑

*n*

_{i}_{=0}−1

*K*(

*x−X*)

_{h}_{n}i*where X*1*= Xt*1*, . . . , Xn= Xtn* *are observed values from a solution to (3.6). The function K*

*is what is called a kernel and hn*a small number named bandwidth. A well-known kernel

*is K(x) = e−x*2/2/√2π*. Note that K*(*x−Xi*

*hn* *)/hn* has area one and is highly concentrated

*around Xifor small hn. An appropriate choice of hnis n*−1/5times the estimated standard

*deviation of X*1*, . . . Xn*. This is somewhat fishy because the standard deviation changes

after the time-change-point.

Nevertheless, the algorithm to simulate change point,θ1 andθ2is as in the case with

known drift:
ˆ
*Zi*=
*Xi*+1*− Xi*
√_{∆}
*n* *− ˆb*
*n(Xi*)
p
∆*n*

*The algorithm is thus as follows: First estimate b by ˆbn*, then apply the algorithm as in

*the case with known drift b replaced by bn*.

There is, as already observed, some things to think about what concerns the case with
*unknown drift. The most serious one is that usually non-parametric estimations ˆbnof b are*

*justified if the the process is stationary, i.e. the distribution of Xt* is the same for all time

*points t. In our case the process is not stationary due to the change in parameters. The*
question of how the algorithm can be theoretically justified in our situation is a concern.

If_{{X}t*: t*≥ 0} is a solution to

*dX* *= b(X)dt +*σ*(X)dW*

*and that solution is stationary, i.e. X*1*= Xt*1*, . . . , Xn= Xtn* have the same density function

π*, one way to show the validity of the approximation ˆbn(x) to b(x) is as follows. First we*

note thatπ*n(x) =*∑*n*_{1}*K((Xi− x)/hn)/(nhn*) is an approximation toπ*(x):*

so
*ˆbn(x) =*
1
*n*∑
*n*
1*h*1*nK*(
*Xi−x*
*hn* *)(Xi*+1*− Xi*)
1
*n*∑
*n*
1*h*1*nK*(
*Xi−x*
*hn* )∆*n*
≈
1
*n*∑
*n*
1*h*1*nK*(
*Xi−x*
*hn* *)(b(Xi*)∆*n*+σ*(Xi*)∆*Wi*)
1
*n*∑
*n*
1
1
*hnK*(
*Xi−x*
*hn* )∆*n*
≈
1
*n*∑
*n*
1*h*1*nK*(
*Xi−x*
*hn* *)b(Xi*)
1
*n*∑
*n*
1*h*1*nK*(
*Xi−x*
*hn* )
≈ *E*[
1
*hnK*(
*X*1*−x*
*hn* *)b(X*1)]
*E*[* _{h}*1

*∑*

_{n}*n*

_{1}

*1*

_{h}*(*

_{n}K*X*1

*−x*

*hn*)] = R

_{1}

*hnK*(

*y−x*

*hn*

*)b(y)*π

*(y)dy*R

_{1}

*hnK*(

*y−x*

*hn*

*)dy*≈

*b(x)*

_{π}

*π*

_{(x)}*(x)*

*= b(x)!*

If_{{X}t*: t* ≥ 0} is as in our case with unknown parameter-shift, it seems not trivial to

*justify the estimate ˆbn* *of b. Nevertheless, the drift estimate also in our situation with*

parameter shift can be justified, [2], [3].

Another question is how to estimate a parameter where the change is not in the
*disper-sion but in the drift, or the change is both in the drift and in the disperdisper-sion as for Zt* given

by (3.7). That is beyond this bachelor thesis. 3.2.1 DWJ

We analyze the DWJ dataset,which contains the weekly closings of the Dow-Jones
in-dustrial average in the period July 1971–August 1974. For references see [2]. There are
162 data and the main evidence found by several authors is that a change in the varience
*occured around k*= 88, which corresponds to the third week of March 19731_{. Instead of}

working on the values we transform the data into returns

*X(ti*) =

*(Y (ti) −Y (ti*−1))

*Y(ti*−1)

*, i = 1, 2, 3...*

*with Y the series of Dow-Jones closings and X the returns. We assume the drift coefficient*
to be unknown. That means we use the method with unknown drift.

See Figures 3.12–3.18. The 1973 change is related to the time when the Watergate scandal became known. A more detailed analysis, publsihed elesewhere, indicates an earlier time change point when there was announced a broken gold/dollar relation. 3.2.2 Goldmann-Sachs data

In the second half of 2008 a world wide finacial crisis started. We wanted to see if we can detect more precisely when that crisis started. Therefore daily closing values from one of the major stock dealers was collected. We chose one stock named ’Goldmann Sachs sr A’ obtained from the database EcoWin, April 25 2005 – May 8, 2009. The computations using the algorithm with unknown drift is used. The figures 3.19–3.25 show a strong indication of time change point. That time point was September 3, 2008.

1_{In March 1973, James McCord wrote a letter to Judge John J. Sirica charging a cover up of the burglary.}

0 0.2 0.4 0.6 0.8 1 −3 −2 −1 0 1 2 3 4 5

*Figure 3.9: Simulation of the Ornstein-Uhlenbeck process X given by the SDE dXt* =

*(6 − 2Xt)dt +*θ*dWt*, whereθ *= 6 for t <*τ0= 0.6 andθ2*= 2 for t ≥ 0.6.X OU*

−3 −2 −1 0 1 2 3 4 5 −40 −30 −20 −10 0 10 20 30 40

*Figure 3.10: The drift b _{(x) = 6 − 2x and its estimate b}n(x) for the Ornstein-Uhlenbeck*

*process X given by the SDE dXt= (6 − 2Xt)dt +*θ*dWt*, whereθ *= 6 for t <*τ0= 0.6 and

0 0.2 0.4 0.6 0.8 1 0 1 2 3 4 5 6 7

Figure 3.11: Estimate ofθ *for the Ornstein-Uhlenbeck process X given by the SDE dXt*=

*(6−2Xt)dt +*θ*dWt*, whereθ *= 6 for t <*τ0= 0.6 andθ2*= 2 for t ≥ 0.6. Estimated change*

time point ˆτ0= 0.5980, ˆθ1= 5.6727, ˆθ2= 1.8679. The estimates ˆθ1 and ˆθ2 are drawn

with a red (dashed) line. The exactθ1andθ2are drawn with a green (solid) line.

1971.5 1972 1972.5 1973 1973.5 1974 1974.5 1975 750 800 850 900 950 1000 1050

1971.5 1972 1972.5 1973 1973.5 1974 1974.5 1975 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08

Figure 3.13: Dow-Jones closings

−0.080 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 2 4 6 8 10 12 14 16 18 20

−0.08−4 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 −3 −2 −1 0 1 2 3 4

*Figure 3.15: Estimate of Dow-Jones drift b*

1971.5 1972 1972.5 1973 1973.5 1974 1974.5 1975 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5

1971.5 1972 1972.5 1973 1973.5 1974 1974.5 1975 0.165 0.17 0.175 0.18 0.185 0.19 0.195 0.2 0.205 0.21 0.215

*Figure 3.17: The Dow-Jones U*2variable. Note the clear change point in the beginning of
1973.
1971.5 1972 1972.5 1973 1973.5 1974 1974.5 1975
0.08
0.1
0.12
0.14
0.16
0.18
0.2
0.22
0.24
0.26
0.28

Figure 3.18: The Dow-Jonesθ variable. Note the clear change point in the beginning of 1973. Estimated change time point ˆτ0= 0.5980, estimated parameters ˆθ1= 0.1100,

ˆ

θ2= 0.2030. The estimates ˆθ1 and ˆθ2 are drawn with a red (dashed) line. In this case

20055 2005.5 2006 2006.5 2007 2007.5 2008 2008.5 2009 2009.5 10 15 20 25 30

Figure 3.19: Goldmann-Sachs closings

2005 2005.5 2006 2006.5 2007 2007.5 2008 2008.5 2009 2009.5 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2

−0.20 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 5 10 15 20 25 30 35 40 45 50

Figure 3.21: Estimate of Goldmann-Sachs returns density. It describes the daily incre-mentw of the Goldmann-Sachs.

−0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 −100 −80 −60 −40 −20 0 20 40 60 80

2005 2005.5 2006 2006.5 2007 2007.5 2008 2008.5 2009 2009.5 −3 −2 −1 0 1 2 3

*Figure 3.23: The Goldmann-Sachs Z variable. Note the change in 2008.*

2005 2005.5 2006 2006.5 2007 2007.5 2008 2008.5 2009 2009.5 390 400 410 420 430 440 450

20050 2005.5 2006 2006.5 2007 2007.5 2008 2008.5 2009 2009.5 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Figure 3.25: The Goldmann-Sachsθ variable. Note the clear change point. It is at Sep 3, 2008. The figure indicates that the financial crisis that appeared 2008 started that particular day. Estimated parameters ˆθ1= 0.1720, ˆθ2= 0.7920. The estimates ˆθ1and ˆθ2

are drawn with a red (dashed) line. In this case there are no exactθ1andθ2: we just have

a model adapted to real data

**4**

**appendix**

*For X given by*

*dXt* *= b(Xt)dt +*σ*(Xt)dW*

*and f a twice differentiable function, it is well known that _{{ f (X}t) : t ≥ 0} also satisfies an*

SDE and is given by the Itô formula

*d f(Xt) = fx(Xt)dXt*+
1
2*fxx(Xt)(dXt*)
2_{,}
where
*d(Xt*)2=σ2*(Xt)dt.*

**4.1** **Typical simulation of an SDE with change-point**

n=1000;
x=1;
T=1;
dt=T/(n-1);
t=linspace(0,T,n);
dW=sqrt(dt)*randn(n,1);
theta=1*(t<1/2)+10*(t>=1/2);
sigma=@(x)1;
for i=1:n-1
x(i+1)=x(i)+b(x(i))*dt+theta(i)*sigma(x(i))*dW(i);
end
**4.2** **Known drift**

Z=(diff(x)-b(x(1:end-1))*dt)/(sqrt(dt)*sigma(x(1:end-1))); S=cumsum(Z.^2); theta1hat=sqrt(S./(1:n-1)); theta2hat=sqrt((S(n-1)-S)./fliplr(1:n-1)); U2=sum((Z.^2-theta2hat.^2).^2)+cumsum((Z.^2-theta1hat.^2... ).^2-(Z.^2-theta2hat.^2).^2); k0hat=find(U2==min(U2)); tau0hat=t(k0hat) thetahat=theta1hat.*(t(1:n-1)<tau0hat)+theta2hat.*(t(1:n-1)>=tau0hat); thetap=theta1hat(k0hat).*(t(1:n-1)<tau0hat)+theta2hat(k0hat).*(t(1:n-1)>=t [t(k0hat) theta1hat(k0hat) theta2hat(k0hat)] ...

%estimates of time shift, theta1

**4.3** **Unknown drift**
h=std(X)*length(X)^(-1/5);
K=@(x)exp(-x.^2/2)/sqrt(2*pi);
pihat=@(x)mean(K((x-X)/h))/h;
bhat=@(x)sum(K((x-X(1:end-1))/h).*diff(X))/(dt*sum(K((x-X(1:end-1))/h)));
sigmahat=@(x)sqrt(sum(K((x-X(1:end-1))/h).*(diff(X)).^2)/(dt*sum(K((x-X(1:
x=linspace(min(X),max(X));
for i=1:length(x)
pihatx(i)=pihat(x(i));
bhatx(i)=bhat(x(i));
end
Z=(diff(X)-bhatX(1:end-1)*dt)./(sqrt(dt)*sigma(X(1:end-1)));
S=cumsum(Z.^2);
theta2hat=sqrt((S(n-1)-S)./fliplr(1:n-1));
U2=sum((Z.^2-theta2hat.^2).^2)+cumsum((Z.^2-theta1hat.^2...
).^2-(Z.^2-theta2hat.^2).^2);
k0hat=find(U2==min(U2));
tau0hat=t(k0hat)
thetahat=theta1hat.*(t(1:n-1)<tau0hat)+theta2hat.*(t(1:n-1)>=tau0hat);
thetap=theta1hat(k0hat).*(t(1:n-1)<tau0hat)+theta2hat(k0hat).*(t(1:n-1)>=t
[t(k0hat) theta1hat(k0hat) theta2hat(k0hat)] ...

%estimates of time shift, theta1

**References**

[1] P. Kloeden, E. Platen, Numerical simulation of stochastic differential equations, Springer, 1999.

[2] S. M. Iacus, Simulation and inference for Stochastic differential equations, Springer, 2008.

## Växjö

## universitet

**Matematiska och systemtekniska institutionen**
SE-351 95 Växjö