Växjö University

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−x2/2/√2π,

x_{∈ R.}

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

(i) Wt is continuous and W0= 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 Fs, 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_{Wt : 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−Wiis a normal (Gaussian) random number with

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

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 X0and a trajectory{Xt} at gridpoint 0 = t0< t1< . . . < 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−Wiis √∆nN(0,1)-distributed, implies that Ziis θ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 Xias 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:_|Zi| given by (3.4) with data Xias 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 Xias simulated in Figure 3.1.θ = 1 for t <τ0= 0.5

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

k0 would be known, then an estimate ofθ1would be simply

ˆ θ1= s ∑k0 1 Zi2 k0

andθ2would be estimated by

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

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

τ0. How do we do then?

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

Sk= k

∑

i=1

Z_i2. (3.5)

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

ˆ θ1(k) = s ∑n i=1Zi2 k = r 1 kSk,

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

What creature is that? U_k2 is a simultaneous measure of the error between ˆθ1(k) and θ1

and between ˆθ2(k) andθ2. If i≤ k0, then Z_i2≈θ₁2and if i> k0, then Z_i2≈θ₂2.

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

ˆk0= arg min k { k

∑

∑

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

∑

∑

∑

∑

∑

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

of U2, ˆθ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 H0:θ1=θ2, then

r n

2|Dk|

d

→ |W0(τ_)|,

where_{W0(τ_{), 0 ≤}τ_{≤ 1} is a Brownian bridge ending and starting at zero.}

By saying ’under the null-hypothesis H0: θ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: U2 for X given by dXt = Xtdt+θdWt, where θ = 1 for t < k0= 0.5 and

θ = 10 for t ≥τ0. Estimated change point ˆτ0= 0.500. For each time point tk the U_k2

-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θ1for 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θ2for 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 < k0= 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= ˆk0/n satisfies

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

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θ₀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_(y)) which can also write as

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_−Xi hn )(Xi+1− Xi) ∆n∑n_i=0−1K(x−X_hni)

where X1= Xt1, . . . , Xn= Xtn are observed values from a solution to (3.6). The function K

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

is K(x) = e−x2/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 X1, . . . 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_{Xt : t≥ 0} is a solution to

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

and that solution is stationary, i.e. X1= Xt1, . . . , 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₁K((Xi− x)/hn)/(nhn) is an approximation toπ(x):

so ˆbn(x) = 1 n∑ n 1h1nK( Xi−x hn )(Xi+1− Xi) 1 n∑ n 1h1nK( Xi−x hn )∆n ≈ 1 n∑ n 1h1nK( Xi−x hn )(b(Xi)∆n+σ(Xi)∆Wi) 1 n∑ n 1 1 hnK( Xi−x hn )∆n ≈ 1 n∑ n 1h1nK( Xi−x hn )b(Xi) 1 n∑ n 1h1nK( Xi−x hn ) ≈ E[ 1 hnK( X1−x hn )b(X1)] E[_h1_n∑n_1h1_nK(X1−x hn )] = R ₁ hnK( y−x hn )b(y)π(y)dy R ₁ hnK( y−x hn )dy ≈ b(x)_π(x)π(x) = b(x)!

If_{Xt : 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 U2variable. 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

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 2fxx(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ö