Quantcast

Ornstein-Uhlenbeck

classic Classic list List threaded Threaded
19 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Ornstein-Uhlenbeck

Stephen Choularton-3
Hi

Wonder if anyone could point me how I use this method to discover the half life of a mean reverting process. 

I am looking into pair trading and the time it takes for a cointegrated pair to revert to the norm.
--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Sarbo
By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to discretise the equation and then fit the data to it. I don't have the time right now to go into all the details but it's not hard- you can parameterise the process using simple linear regression. If you need help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
Hi

Wonder if anyone could point me how I use this method to discover the half life of a mean reverting process. 

I am looking into pair trading and the time it takes for a cointegrated pair to revert to the norm.
--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Stale-2
In reply to this post by Stephen Choularton-3
Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn






> ------------------------------
>
> Message: 2
> Date: Tue, 12 Oct 2010 05:43:32 -0400
> From: Sarbo <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="utf-8"
>
> By half-life, do you mean the speed of mean-reversion?
>
> If so, there's a bit of algebraic tomfoolery that's required to
> discretise the equation and then fit the data to it. I don't have the
> time right now to go into all the details but it's not hard- you can
> parameterise the process using simple linear regression. If you need
> help with that I'll try and get back to you tonight about it.
>
> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>
> > Hi
> >
> > Wonder if anyone could point me how I use this method to discover the
> > half life of a mean reverting process.
> >
> > I am looking into pair trading and the time it takes for a
> > cointegrated pair to revert to the norm.
> >
> > --
> > Stephen Choularton Ph.D., FIoD
> >
> > 9999 2226
> > 0413 545 182
> >
> >
> > for insurance go to www.netinsure.com.au
> > for markets go to www.organicfoodmarkets.com.au
> >
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R questions
> should go.
>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
> >
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: CoS2010Winner.JPG
> Type: image/jpeg
> Size: 16091 bytes
> Desc: not available
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
> >
>
> ------------------------------
>
> _______________________________________________
> R-SIG-Finance mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>
> End of R-SIG-Finance Digest, Vol 77, Issue 8
> ********************************************
        [[alternative HTML version deleted]]


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Paul Teetor
Stephen:

I do mean-reversion trading, and I use a half-life analysis to judge the
wisdom of a trade. If the estimated half-life is too long, it doesn't make
sense to take the trade. It's a time-vs-risk thing.

In the past, I used the log(2)/speed formula mentioned by Bjorn, below. (His
link is very useful, BTW.) However, I was very unhappy with the estimates
provided by that formula. They did not match my actual trading experience.

I did some research on the topic, and got some useful results. I added a
momentum term to my model, measuruing the current slope of the spread. The
slope answers an important question: is the spread currently reverting
(moving towards the mean) or is it averting (moving away from the mean)? The
half-life is different, depending upon the current phase (reversion vs.
aversion). I found this conditioning term was statistically significant, so
I condition my estimate on it. I generate historical data and partition it
according to the market's reversion/aversion state. The two partitions shows
a different half-life, with the reverting phase having a (much) shorter
half-life than the averting phase. Both phases show an exponentally
distributed half-life, but with different means. I use those historical
estimates now in my trading, instead of the O-U estimate.

In retrospect, the problem with the O-U model is that it assumes the
mean-reverting process is always reverting. That is not quite correct. In my
experience, the spreads can alternate between periods of mean-aversion and
mean-reversion.

I hope that helps (and I hope it makes sense!).

Paul


-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Bjorn Skogtro
Sent: Tuesday, October 12, 2010 5:34 AM
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck

Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjxrn






> ------------------------------
>
> Message: 2
> Date: Tue, 12 Oct 2010 05:43:32 -0400
> From: Sarbo <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="utf-8"
>
> By half-life, do you mean the speed of mean-reversion?
>
> If so, there's a bit of algebraic tomfoolery that's required to
> discretise the equation and then fit the data to it. I don't have the
> time right now to go into all the details but it's not hard- you can
> parameterise the process using simple linear regression. If you need
> help with that I'll try and get back to you tonight about it.
>
> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>
> > Hi
> >
> > Wonder if anyone could point me how I use this method to discover
> > the half life of a mean reverting process.
> >
> > I am looking into pair trading and the time it takes for a
> > cointegrated pair to revert to the norm.
> >
> > --
> > Stephen Choularton Ph.D., FIoD
> >
> > 9999 2226
> > 0413 545 182
> >
> >
> > for insurance go to www.netinsure.com.au for markets go to
> > www.organicfoodmarkets.com.au
> >
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R
> > questions
> should go.
>
>
> -------------- next part -------------- An HTML attachment was
> scrubbed...
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e3
> 2fc7/attachment-0001.html
> >
> -------------- next part -------------- A non-text attachment was
> scrubbed...
> Name: CoS2010Winner.JPG
> Type: image/jpeg
> Size: 16091 bytes
> Desc: not available
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e3
> 2fc7/attachment-0001.jpe
> >
>
> ------------------------------
>
> _______________________________________________
> R-SIG-Finance mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>
> End of R-SIG-Finance Digest, Vol 77, Issue 8
> ********************************************

        [[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

stefano iacus-2
In reply to this post by Stale-2
just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
 }

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)
 
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:

> Hi Stephen,
>
> You could take a look at
>
> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>
> for the linear regression method, or take a look at the package "sde" which
> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
> though, only the CIR).
>
> The half-life is given as log(2)/mean-reversion speed.
>
> Do keep an eye on the partition of the time-axis, e.g. what frequency you
> are using (daily, yearly) for interpreting the half-life.
>
> BR,
> Bjørn
>
>
>
>
>
>
>> ------------------------------
>>
>> Message: 2
>> Date: Tue, 12 Oct 2010 05:43:32 -0400
>> From: Sarbo <[hidden email]>
>> To: [hidden email]
>> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>> Message-ID: <[hidden email]>
>> Content-Type: text/plain; charset="utf-8"
>>
>> By half-life, do you mean the speed of mean-reversion?
>>
>> If so, there's a bit of algebraic tomfoolery that's required to
>> discretise the equation and then fit the data to it. I don't have the
>> time right now to go into all the details but it's not hard- you can
>> parameterise the process using simple linear regression. If you need
>> help with that I'll try and get back to you tonight about it.
>>
>> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>>
>>> Hi
>>>
>>> Wonder if anyone could point me how I use this method to discover the
>>> half life of a mean reverting process.
>>>
>>> I am looking into pair trading and the time it takes for a
>>> cointegrated pair to revert to the norm.
>>>
>>> --
>>> Stephen Choularton Ph.D., FIoD
>>>
>>> 9999 2226
>>> 0413 545 182
>>>
>>>
>>> for insurance go to www.netinsure.com.au
>>> for markets go to www.organicfoodmarkets.com.au
>>>
>>>
>>> _______________________________________________
>>> [hidden email] mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>> -- Subscriber-posting only. If you want to post, subscribe first.
>>> -- Also note that this is not the r-help list where general R questions
>> should go.
>>
>>
>> -------------- next part --------------
>> An HTML attachment was scrubbed...
>> URL: <
>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>>>
>> -------------- next part --------------
>> A non-text attachment was scrubbed...
>> Name: CoS2010Winner.JPG
>> Type: image/jpeg
>> Size: 16091 bytes
>> Desc: not available
>> URL: <
>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>>>
>>
>> ------------------------------
>>
>> _______________________________________________
>> R-SIG-Finance mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>
>>
>> End of R-SIG-Finance Digest, Vol 77, Issue 8
>> ********************************************
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.


-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Stale-2
In reply to this post by Paul Teetor
Just another thought:

If you google around, you may stumble upon some results on first passage
times for the OU-process that may be useful to you. These are known in an
closed-form solution.

I think, like Paul, that polishing the strategy is essential. Finding
reverting pairs are the easy part, but optimizing the entry/exit signals are
the trick.

Bjørn

PS: A simple code for calibrating the OU-process. Should've attached this
the first time :)

#source("ouFit.R")
ouFit=function(spread) {
  n=length(spread)
  x=spread[1:(n-1)]
  y=spread[2:n]
  spread.fit=lm(y~x)
  coefs=as.numeric(coefficients(spread.fit))
  a=coefs[1]
  b=coefs[2]
  err=var(residuals(spread.fit))

  alpha=-log(a)
  mu=alpha*b/(1-a)
  sigma=sqrt(2*alpha*err/(1-a^2))
  theta=list(alpha=alpha, mu=mu, sigma=sigma)
  return(theta)
}

It won't give you half-life, but you can use the alpha from the above code
to find it, in combination with my previous post.


2010/10/12 Paul Teetor <[hidden email]>

> Stephen:
>
> I do mean-reversion trading, and I use a half-life analysis to judge the
> wisdom of a trade. If the estimated half-life is too long, it doesn't make
> sense to take the trade. It's a time-vs-risk thing.
>
> In the past, I used the log(2)/speed formula mentioned by Bjorn, below.
> (His
> link is very useful, BTW.) However, I was very unhappy with the estimates
> provided by that formula. They did not match my actual trading experience.
>
> I did some research on the topic, and got some useful results. I added a
> momentum term to my model, measuruing the current slope of the spread. The
> slope answers an important question: is the spread currently reverting
> (moving towards the mean) or is it averting (moving away from the mean)?
> The
> half-life is different, depending upon the current phase (reversion vs.
> aversion). I found this conditioning term was statistically significant, so
> I condition my estimate on it. I generate historical data and partition it
> according to the market's reversion/aversion state. The two partitions
> shows
> a different half-life, with the reverting phase having a (much) shorter
> half-life than the averting phase. Both phases show an exponentally
> distributed half-life, but with different means. I use those historical
> estimates now in my trading, instead of the O-U estimate.
>
> In retrospect, the problem with the O-U model is that it assumes the
> mean-reverting process is always reverting. That is not quite correct. In
> my
> experience, the spreads can alternate between periods of mean-aversion and
> mean-reversion.
>
> I hope that helps (and I hope it makes sense!).
>
> Paul
>
>
> -----Original Message-----
> From: [hidden email]
> [mailto:[hidden email]] On Behalf Of Bjorn
> Skogtro
> Sent: Tuesday, October 12, 2010 5:34 AM
> To: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>
> Hi Stephen,
>
> You could take a look at
>
> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>
> for the linear regression method, or take a look at the package "sde" which
> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
> though, only the CIR).
>
> The half-life is given as log(2)/mean-reversion speed.
>
> Do keep an eye on the partition of the time-axis, e.g. what frequency you
> are using (daily, yearly) for interpreting the half-life.
>
> BR,
> Bjxrn
>
>
>
>
>
>
> > ------------------------------
> >
> > Message: 2
> > Date: Tue, 12 Oct 2010 05:43:32 -0400
> > From: Sarbo <[hidden email]>
> > To: [hidden email]
> > Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> > Message-ID: <[hidden email]>
> > Content-Type: text/plain; charset="utf-8"
> >
> > By half-life, do you mean the speed of mean-reversion?
> >
> > If so, there's a bit of algebraic tomfoolery that's required to
> > discretise the equation and then fit the data to it. I don't have the
> > time right now to go into all the details but it's not hard- you can
> > parameterise the process using simple linear regression. If you need
> > help with that I'll try and get back to you tonight about it.
> >
> > On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
> >
> > > Hi
> > >
> > > Wonder if anyone could point me how I use this method to discover
> > > the half life of a mean reverting process.
> > >
> > > I am looking into pair trading and the time it takes for a
> > > cointegrated pair to revert to the norm.
> > >
> > > --
> > > Stephen Choularton Ph.D., FIoD
> > >
> > > 9999 2226
> > > 0413 545 182
> > >
> > >
> > > for insurance go to www.netinsure.com.au for markets go to
> > > www.organicfoodmarkets.com.au
> > >
> > >
> > > _______________________________________________
> > > [hidden email] mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > > -- Subscriber-posting only. If you want to post, subscribe first.
> > > -- Also note that this is not the r-help list where general R
> > > questions
> > should go.
> >
> >
> > -------------- next part -------------- An HTML attachment was
> > scrubbed...
> > URL: <
> > https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e3
> > 2fc7/attachment-0001.html
> > >
> > -------------- next part -------------- A non-text attachment was
> > scrubbed...
> > Name: CoS2010Winner.JPG
> > Type: image/jpeg
> > Size: 16091 bytes
> > Desc: not available
> > URL: <
> > https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e3
> > 2fc7/attachment-0001.jpe
> > >
> >
> > ------------------------------
> >
> > _______________________________________________
> > R-SIG-Finance mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> >
> >
> > End of R-SIG-Finance Digest, Vol 77, Issue 8
> > ********************************************
>
>         [[alternative HTML version deleted]]
>
>
>

--
Up, down, turn around
Please dont let me hit the ground
Tonight I think Ill walk alone
Ill find my soul as I go home.

        [[alternative HTML version deleted]]


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Patrick Burns-2
In reply to this post by stefano iacus-2
The OU process is Gaussian, but the
market didn't get the memo that *it*
has to be Gaussian.

On 12/10/2010 16:41, stefano iacus wrote:

> just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
>
> sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
>
> This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
>
>
> # ex3.01.R
> OU.lik<- function(theta1, theta2, theta3){
>    n<- length(X)
>    dt<- deltat(X)
>    -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
>   }
>
> require(stats4)
> require(sde)
> set.seed(123)
> X<- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>        method="L-BFGS-B", lower=c(-Inf,0,0)) ->  fit
> summary(fit)
>
> # ex3.01.R (cont.)
> prof<- profile(fit)
> par(mfrow=c(1,3))
> plot(prof)
> par(mfrow=c(1,1))
> vcov(fit)
> confint(fit)
>
> # ex3.01.R (cont.)
> set.seed(123)
> X<- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>        method="L-BFGS-B", lower=c(-Inf,0,0)) ->  fit2
> summary(fit2)
>
>
> I hope this helps out
>
> stefano
>
> On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
>
>> Hi Stephen,
>>
>> You could take a look at
>>
>> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>>
>> for the linear regression method, or take a look at the package "sde" which
>> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
>> though, only the CIR).
>>
>> The half-life is given as log(2)/mean-reversion speed.
>>
>> Do keep an eye on the partition of the time-axis, e.g. what frequency you
>> are using (daily, yearly) for interpreting the half-life.
>>
>> BR,
>> Bjørn
>>
>>
>>
>>
>>
>>
>>> ------------------------------
>>>
>>> Message: 2
>>> Date: Tue, 12 Oct 2010 05:43:32 -0400
>>> From: Sarbo<[hidden email]>
>>> To: [hidden email]
>>> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>>> Message-ID:<[hidden email]>
>>> Content-Type: text/plain; charset="utf-8"
>>>
>>> By half-life, do you mean the speed of mean-reversion?
>>>
>>> If so, there's a bit of algebraic tomfoolery that's required to
>>> discretise the equation and then fit the data to it. I don't have the
>>> time right now to go into all the details but it's not hard- you can
>>> parameterise the process using simple linear regression. If you need
>>> help with that I'll try and get back to you tonight about it.
>>>
>>> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>>>
>>>> Hi
>>>>
>>>> Wonder if anyone could point me how I use this method to discover the
>>>> half life of a mean reverting process.
>>>>
>>>> I am looking into pair trading and the time it takes for a
>>>> cointegrated pair to revert to the norm.
>>>>
>>>> --
>>>> Stephen Choularton Ph.D., FIoD
>>>>
>>>> 9999 2226
>>>> 0413 545 182
>>>>
>>>>
>>>> for insurance go to www.netinsure.com.au
>>>> for markets go to www.organicfoodmarkets.com.au
>>>>
>>>>
>>>> _______________________________________________
>>>> [hidden email] mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>>> -- Subscriber-posting only. If you want to post, subscribe first.
>>>> -- Also note that this is not the r-help list where general R questions
>>> should go.
>>>
>>>
>>> -------------- next part --------------
>>> An HTML attachment was scrubbed...
>>> URL:<
>>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>>>>
>>> -------------- next part --------------
>>> A non-text attachment was scrubbed...
>>> Name: CoS2010Winner.JPG
>>> Type: image/jpeg
>>> Size: 16091 bytes
>>> Desc: not available
>>> URL:<
>>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>>>>
>>>
>>> ------------------------------
>>>
>>> _______________________________________________
>>> R-SIG-Finance mailing list
>>> [hidden email]
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>>
>>>
>>> End of R-SIG-Finance Digest, Vol 77, Issue 8
>>> ********************************************
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions should go.
>
>
> -----------------------------------
> Stefano M. Iacus
> Department of Economics,
> Business and Statistics
> University of Milan
> Via Conservatorio, 7
> I-20123 Milan - Italy
> Ph.: +39 02 50321 461
> Fax: +39 02 50321 505
> http://www.economia.unimi.it/iacus
> ------------------------------------------------------------------------------------
> Please don't send me Word or PowerPoint attachments if not
> absolutely necessary. See:
> http://www.gnu.org/philosophy/no-word-attachments.html
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>

--
Patrick Burns
[hidden email]
http://www.burns-stat.com
http://www.portfolioprobe.com/blog

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

stefano iacus-2
that's another point and I completely agree with you :-)

so the real point is not trying to fit wrong models on the (so to say) "wrong" data


stefano

On 12 Oct 2010, at 18:00, Patrick Burns wrote:

> The OU process is Gaussian, but the
> market didn't get the memo that *it*
> has to be Gaussian.
>
> On 12/10/2010 16:41, stefano iacus wrote:
>> just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
>>
>> sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
>>
>> This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
>>
>>
>> # ex3.01.R
>> OU.lik<- function(theta1, theta2, theta3){
>>   n<- length(X)
>>   dt<- deltat(X)
>>   -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
>>  }
>>
>> require(stats4)
>> require(sde)
>> set.seed(123)
>> X<- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>>       method="L-BFGS-B", lower=c(-Inf,0,0)) ->  fit
>> summary(fit)
>>
>> # ex3.01.R (cont.)
>> prof<- profile(fit)
>> par(mfrow=c(1,3))
>> plot(prof)
>> par(mfrow=c(1,1))
>> vcov(fit)
>> confint(fit)
>>
>> # ex3.01.R (cont.)
>> set.seed(123)
>> X<- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>>       method="L-BFGS-B", lower=c(-Inf,0,0)) ->  fit2
>> summary(fit2)
>>
>>
>> I hope this helps out
>>
>> stefano
>>
>> On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
>>
>>> Hi Stephen,
>>>
>>> You could take a look at
>>>
>>> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>>>
>>> for the linear regression method, or take a look at the package "sde" which
>>> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
>>> though, only the CIR).
>>>
>>> The half-life is given as log(2)/mean-reversion speed.
>>>
>>> Do keep an eye on the partition of the time-axis, e.g. what frequency you
>>> are using (daily, yearly) for interpreting the half-life.
>>>
>>> BR,
>>> Bjørn
>>>
>>>
>>>
>>>
>>>
>>>
>>>> ------------------------------
>>>>
>>>> Message: 2
>>>> Date: Tue, 12 Oct 2010 05:43:32 -0400
>>>> From: Sarbo<[hidden email]>
>>>> To: [hidden email]
>>>> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>>>> Message-ID:<[hidden email]>
>>>> Content-Type: text/plain; charset="utf-8"
>>>>
>>>> By half-life, do you mean the speed of mean-reversion?
>>>>
>>>> If so, there's a bit of algebraic tomfoolery that's required to
>>>> discretise the equation and then fit the data to it. I don't have the
>>>> time right now to go into all the details but it's not hard- you can
>>>> parameterise the process using simple linear regression. If you need
>>>> help with that I'll try and get back to you tonight about it.
>>>>
>>>> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>>>>
>>>>> Hi
>>>>>
>>>>> Wonder if anyone could point me how I use this method to discover the
>>>>> half life of a mean reverting process.
>>>>>
>>>>> I am looking into pair trading and the time it takes for a
>>>>> cointegrated pair to revert to the norm.
>>>>>
>>>>> --
>>>>> Stephen Choularton Ph.D., FIoD
>>>>>
>>>>> 9999 2226
>>>>> 0413 545 182
>>>>>
>>>>>
>>>>> for insurance go to www.netinsure.com.au
>>>>> for markets go to www.organicfoodmarkets.com.au
>>>>>
>>>>>
>>>>> _______________________________________________
>>>>> [hidden email] mailing list
>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>>>> -- Subscriber-posting only. If you want to post, subscribe first.
>>>>> -- Also note that this is not the r-help list where general R questions
>>>> should go.
>>>>
>>>>
>>>> -------------- next part --------------
>>>> An HTML attachment was scrubbed...
>>>> URL:<
>>>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>>>>>
>>>> -------------- next part --------------
>>>> A non-text attachment was scrubbed...
>>>> Name: CoS2010Winner.JPG
>>>> Type: image/jpeg
>>>> Size: 16091 bytes
>>>> Desc: not available
>>>> URL:<
>>>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>>>>>
>>>>
>>>> ------------------------------
>>>>
>>>> _______________________________________________
>>>> R-SIG-Finance mailing list
>>>> [hidden email]
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>>>
>>>>
>>>> End of R-SIG-Finance Digest, Vol 77, Issue 8
>>>> ********************************************
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> [hidden email] mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>> -- Subscriber-posting only. If you want to post, subscribe first.
>>> -- Also note that this is not the r-help list where general R questions should go.
>>
>>
>> -----------------------------------
>> Stefano M. Iacus
>> Department of Economics,
>> Business and Statistics
>> University of Milan
>> Via Conservatorio, 7
>> I-20123 Milan - Italy
>> Ph.: +39 02 50321 461
>> Fax: +39 02 50321 505
>> http://www.economia.unimi.it/iacus
>> ------------------------------------------------------------------------------------
>> Please don't send me Word or PowerPoint attachments if not
>> absolutely necessary. See:
>> http://www.gnu.org/philosophy/no-word-attachments.html
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions should go.
>>
>
> --
> Patrick Burns
> [hidden email]
> http://www.burns-stat.com
> http://www.portfolioprobe.com/blog
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.


-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Stephen Choularton-3
In reply to this post by stefano iacus-2
Thanks for this help.

Trying to make sense of it so I have added some notes to the code.  I have marked them #?#

Delighted if you can tell me if I am write or wrong, add any comments, answers.

#?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck 
#?# process work' particularly via dcOU
#?# I have noted in several places that I am after:
#?# 	'the half-life of the decay equals ln(2)/θ'
#?# 	'The half-life is given as log(2)/mean-reversion speed.'
#?#  and I see theta appearing at a number of points in the code.
#?#  Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
#?#  eg is one of these the theta I am after?

# ex3.01.R 
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
 }

require(stats4)
require(sde)

#?# random numer generation seed
set.seed(123)

#?# creation of a data set
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
#?# If I Look at X its like this:
#?# Time Series:
#?# Start = 0 
#?# End = 1000 
#?# Frequency = 1 
#?# [1]  1.00000000  etc
#?# What sort of data object is it and how would I coerce an object with one 
#?# column from a read.csv into it? 


mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

#?# This gives:

#?# Maximum likelihood estimation

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5, 
#?#     theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?#        Estimate Std. Error
#?# theta1 3.355322 0.28159504
#?# theta2 1.106107 0.09010627
#?# theta3 2.052815 0.07624441

#?# -2 log L: 3366.389 

#?# What's this telling me?

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

#?# This provides me with this output using 'fit' from before:

#?# > vcov(fit)
#?#            theta1      theta2      theta3
#?# theta1 0.07929576 0.024620718 0.016634557
#?# theta2 0.02462072 0.008119141 0.005485549
#?# theta3 0.01663456 0.005485549 0.005813209
#?# > confint(fit)
#?# Profiling...
#?#            2.5 %   97.5 %
#?# theta1 2.8448980 3.960982
#?# theta2 0.9433338 1.300629
#?# theta3 1.9147136 2.216113

#?# and 'fit' is:

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5, 
#?#     theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?#   theta1   theta2   theta3 
#?# 3.355322 1.106107 2.052815 

#?# plus some graphic output

#?# Again, what's this telling me.

#?# This looks like a further example?
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)




Please excuse the length of this email (and my lack of understanding)

Hope you can help and thanks.



Stephen Choularton Ph.D., FIoD


On 13/10/2010 2:41 AM, stefano iacus wrote:
just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
 }

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)
 
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:

Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn






------------------------------

Message: 2
Date: Tue, 12 Oct 2010 05:43:32 -0400
From: Sarbo [hidden email]
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
Message-ID: [hidden email]
Content-Type: text/plain; charset="utf-8"

By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to
discretise the equation and then fit the data to it. I don't have the
time right now to go into all the details but it's not hard- you can
parameterise the process using simple linear regression. If you need
help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:

Hi

Wonder if anyone could point me how I use this method to discover the
half life of a mean reverting process.

I am looking into pair trading and the time it takes for a
cointegrated pair to revert to the norm.

--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions
should go.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html

          
-------------- next part --------------
A non-text attachment was scrubbed...
Name: CoS2010Winner.JPG
Type: image/jpeg
Size: 16091 bytes
Desc: not available
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe

          
------------------------------

_______________________________________________
R-SIG-Finance mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-finance


End of R-SIG-Finance Digest, Vol 77, Issue 8
********************************************
	[[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.

-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not 
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.

No virus found in this incoming message. Checked by AVG - www.avg.com

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Ornstein-Uhlenbeck

Stephen Choularton-3
Hi

I am still trying to sort this one out.  Any comments from anyone would be most welcome.

Stephen Choularton Ph.D., FIoD



On 14/10/2010 7:29 AM, Stephen Choularton wrote:
Thanks for this help.

Trying to make sense of it so I have added some notes to the code.  I have marked them #?#

Delighted if you can tell me if I am write or wrong, add any comments, answers.

#?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck 
#?# process work' particularly via dcOU
#?# I have noted in several places that I am after:
#?# 	'the half-life of the decay equals ln(2)/θ'
#?# 	'The half-life is given as log(2)/mean-reversion speed.'
#?#  and I see theta appearing at a number of points in the code.
#?#  Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
#?#  eg is one of these the theta I am after?

# ex3.01.R 
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
 }

require(stats4)
require(sde)

#?# random numer generation seed
set.seed(123)

#?# creation of a data set
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
#?# If I Look at X its like this:
#?# Time Series:
#?# Start = 0 
#?# End = 1000 
#?# Frequency = 1 
#?# [1]  1.00000000  etc
#?# What sort of data object is it and how would I coerce an object with one 
#?# column from a read.csv into it? 


mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

#?# This gives:

#?# Maximum likelihood estimation

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5, 
#?#     theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?#        Estimate Std. Error
#?# theta1 3.355322 0.28159504
#?# theta2 1.106107 0.09010627
#?# theta3 2.052815 0.07624441

#?# -2 log L: 3366.389 

#?# What's this telling me?

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

#?# This provides me with this output using 'fit' from before:

#?# > vcov(fit)
#?#            theta1      theta2      theta3
#?# theta1 0.07929576 0.024620718 0.016634557
#?# theta2 0.02462072 0.008119141 0.005485549
#?# theta3 0.01663456 0.005485549 0.005813209
#?# > confint(fit)
#?# Profiling...
#?#            2.5 %   97.5 %
#?# theta1 2.8448980 3.960982
#?# theta2 0.9433338 1.300629
#?# theta3 1.9147136 2.216113

#?# and 'fit' is:

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5, 
#?#     theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?#   theta1   theta2   theta3 
#?# 3.355322 1.106107 2.052815 

#?# plus some graphic output

#?# Again, what's this telling me.

#?# This looks like a further example?
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)




Please excuse the length of this email (and my lack of understanding)

Hope you can help and thanks.



Stephen Choularton Ph.D., FIoD


On 13/10/2010 2:41 AM, stefano iacus wrote:
just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
  n <- length(X)
  dt <- deltat(X)
  -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
 }

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)
 
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1), 
      method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:

Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn






------------------------------

Message: 2
Date: Tue, 12 Oct 2010 05:43:32 -0400
From: Sarbo [hidden email]
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
Message-ID: [hidden email]
Content-Type: text/plain; charset="utf-8"

By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to
discretise the equation and then fit the data to it. I don't have the
time right now to go into all the details but it's not hard- you can
parameterise the process using simple linear regression. If you need
help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:

Hi

Wonder if anyone could point me how I use this method to discover the
half life of a mean reverting process.

I am looking into pair trading and the time it takes for a
cointegrated pair to revert to the norm.

--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions
should go.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: CoS2010Winner.JPG
Type: image/jpeg
Size: 16091 bytes
Desc: not available
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
------------------------------

_______________________________________________
R-SIG-Finance mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-finance


End of R-SIG-Finance Digest, Vol 77, Issue 8
********************************************
	[[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not 
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.

No virus found in this incoming message. Checked by AVG - www.avg.com
_______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-finance -- Subscriber-posting only. If you want to post, subscribe first. -- Also note that this is not the r-help list where general R questions should go.
No virus found in this incoming message. Checked by AVG - www.avg.com

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Mean reversion

Yihao Lu aeolus_lu

I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
Is there anyone who can give me some possible explanation or guide me to some reference? thanks

Best,
Yihao







________________________________

> Date: Tue, 19 Oct 2010 09:03:55 +1100
> From: [hidden email]
> To: [hidden email]
> CC: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>
> Hi
>
> I am still trying to sort this one out. Any comments from anyone would
> be most welcome.
>
> Stephen Choularton Ph.D., FIoD
>
>
>
> On 14/10/2010 7:29 AM, Stephen Choularton wrote:
> Thanks for this help.
>
> Trying to make sense of it so I have added some notes to the code. I
> have marked them #?#
>
> Delighted if you can tell me if I am write or wrong, add any comments,
> answers.
>
>
> #?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
> #?# process work' particularly via dcOU
> #?# I have noted in several places that I am after:
> #?# 'the half-life of the decay equals ln(2)/θ'
> #?# 'The half-life is given as log(2)/mean-reversion speed.'
> #?# and I see theta appearing at a number of points in the code.
> #?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
> #?# eg is one of these the theta I am after?
>
> # ex3.01.R
> OU.lik <- function(theta1, theta2, theta3){
> n <- length(X)
> dt <- deltat(X)
> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> }
>
> require(stats4)
> require(sde)
>
> #?# random numer generation seed
> set.seed(123)
>
> #?# creation of a data set
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> #?# If I Look at X its like this:
> #?# Time Series:
> #?# Start = 0
> #?# End = 1000
> #?# Frequency = 1
> #?# [1] 1.00000000 etc
> #?# What sort of data object is it and how would I coerce an object with one
> #?# column from a read.csv into it?
>
>
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> summary(fit)
>
> #?# This gives:
>
> #?# Maximum likelihood estimation
>
> #?# Call:
> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>
> #?# Coefficients:
> #?# Estimate Std. Error
> #?# theta1 3.355322 0.28159504
> #?# theta2 1.106107 0.09010627
> #?# theta3 2.052815 0.07624441
>
> #?# -2 log L: 3366.389
>
> #?# What's this telling me?
>
> # ex3.01.R (cont.)
> prof <- profile(fit)
> par(mfrow=c(1,3))
> plot(prof)
> par(mfrow=c(1,1))
> vcov(fit)
> confint(fit)
>
> #?# This provides me with this output using 'fit' from before:
>
> #?# > vcov(fit)
> #?# theta1 theta2 theta3
> #?# theta1 0.07929576 0.024620718 0.016634557
> #?# theta2 0.02462072 0.008119141 0.005485549
> #?# theta3 0.01663456 0.005485549 0.005813209
> #?# > confint(fit)
> #?# Profiling...
> #?# 2.5 % 97.5 %
> #?# theta1 2.8448980 3.960982
> #?# theta2 0.9433338 1.300629
> #?# theta3 1.9147136 2.216113
>
> #?# and 'fit' is:
>
> #?# Call:
> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>
> #?# Coefficients:
> #?# theta1 theta2 theta3
> #?# 3.355322 1.106107 2.052815
>
> #?# plus some graphic output
>
> #?# Again, what's this telling me.
>
> #?# This looks like a further example?
> # ex3.01.R (cont.)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> summary(fit2)
>
>
>
>
> Please excuse the length of this email (and my lack of understanding)
>
> Hope you can help and thanks.
>
>
>
>
> Stephen Choularton Ph.D., FIoD
>
>
> On 13/10/2010 2:41 AM, stefano iacus wrote:
>
> just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
>
> sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
>
> This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
>
>
> # ex3.01.R
> OU.lik <- function(theta1, theta2, theta3){
> n <- length(X)
> dt <- deltat(X)
> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> }
>
> require(stats4)
> require(sde)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> summary(fit)
>
> # ex3.01.R (cont.)
> prof <- profile(fit)
> par(mfrow=c(1,3))
> plot(prof)
> par(mfrow=c(1,1))
> vcov(fit)
> confint(fit)
>
> # ex3.01.R (cont.)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> summary(fit2)
>
>
> I hope this helps out
>
> stefano
>
> On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
>
>
>
> Hi Stephen,
>
> You could take a look at
>
> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>
> for the linear regression method, or take a look at the package "sde" which
> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
> though, only the CIR).
>
> The half-life is given as log(2)/mean-reversion speed.
>
> Do keep an eye on the partition of the time-axis, e.g. what frequency you
> are using (daily, yearly) for interpreting the half-life.
>
> BR,
> Bjørn
>
>
>
>
>
>
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 12 Oct 2010 05:43:32 -0400
> From: Sarbo
> To: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> Message-ID:
> Content-Type: text/plain; charset="utf-8"
>
> By half-life, do you mean the speed of mean-reversion?
>
> If so, there's a bit of algebraic tomfoolery that's required to
> discretise the equation and then fit the data to it. I don't have the
> time right now to go into all the details but it's not hard- you can
> parameterise the process using simple linear regression. If you need
> help with that I'll try and get back to you tonight about it.
>
> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>
>
>
> Hi
>
> Wonder if anyone could point me how I use this method to discover the
> half life of a mean reverting process.
>
> I am looking into pair trading and the time it takes for a
> cointegrated pair to revert to the norm.
>
> --
> Stephen Choularton Ph.D., FIoD
>
> 9999 2226
> 0413 545 182
>
>
> for insurance go to www.netinsure.com.au
> for markets go to www.organicfoodmarkets.com.au
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions
>
>
> should go.
>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>
>
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: CoS2010Winner.JPG
> Type: image/jpeg
> Size: 16091 bytes
> Desc: not available
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>
>
> ------------------------------
>
> _______________________________________________
> R-SIG-Finance mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>
> End of R-SIG-Finance Digest, Vol 77, Issue 8
> ********************************************
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
> -----------------------------------
> Stefano M. Iacus
> Department of Economics,
> Business and Statistics
> University of Milan
> Via Conservatorio, 7
> I-20123 Milan - Italy
> Ph.: +39 02 50321 461
> Fax: +39 02 50321 505
> http://www.economia.unimi.it/iacus
> ------------------------------------------------------------------------------------
> Please don't send me Word or PowerPoint attachments if not
> absolutely necessary. See:
> http://www.gnu.org/philosophy/no-word-attachments.html
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
>
>
>
>
>
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
>
>
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
> Subscriber-posting only. If you want to post, subscribe first. -- Also
> note that this is not the r-help list where general R questions should
> go.
     
_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Mean reversion

Stephen Choularton-3
Wish I could, but as you can see I am having even greater problems ;-)

Stephen Choularton Ph.D., FIoD

On 19/10/2010 12:35 PM, Yihao Lu aeolus_lu wrote:
I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
Is there anyone who can give me some possible explanation or guide me to some reference? thanks

Best,
Yihao







________________________________
Date: Tue, 19 Oct 2010 09:03:55 +1100
From: [hidden email]
To: [hidden email]
CC: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck

Hi

I am still trying to sort this one out. Any comments from anyone would
be most welcome.

Stephen Choularton Ph.D., FIoD



On 14/10/2010 7:29 AM, Stephen Choularton wrote:
Thanks for this help.

Trying to make sense of it so I have added some notes to the code. I
have marked them #?#

Delighted if you can tell me if I am write or wrong, add any comments,
answers.


#?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
#?# process work' particularly via dcOU
#?# I have noted in several places that I am after:
#?# 'the half-life of the decay equals ln(2)/θ'
#?# 'The half-life is given as log(2)/mean-reversion speed.'
#?# and I see theta appearing at a number of points in the code.
#?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
#?# eg is one of these the theta I am after?

# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)

#?# random numer generation seed
set.seed(123)

#?# creation of a data set
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
#?# If I Look at X its like this:
#?# Time Series:
#?# Start = 0
#?# End = 1000
#?# Frequency = 1
#?# [1] 1.00000000 etc
#?# What sort of data object is it and how would I coerce an object with one
#?# column from a read.csv into it?


mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

#?# This gives:

#?# Maximum likelihood estimation

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# Estimate Std. Error
#?# theta1 3.355322 0.28159504
#?# theta2 1.106107 0.09010627
#?# theta3 2.052815 0.07624441

#?# -2 log L: 3366.389

#?# What's this telling me?

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

#?# This provides me with this output using 'fit' from before:

#?# > vcov(fit)
#?# theta1 theta2 theta3
#?# theta1 0.07929576 0.024620718 0.016634557
#?# theta2 0.02462072 0.008119141 0.005485549
#?# theta3 0.01663456 0.005485549 0.005813209
#?# > confint(fit)
#?# Profiling...
#?# 2.5 % 97.5 %
#?# theta1 2.8448980 3.960982
#?# theta2 0.9433338 1.300629
#?# theta3 1.9147136 2.216113

#?# and 'fit' is:

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# theta1 theta2 theta3
#?# 3.355322 1.106107 2.052815

#?# plus some graphic output

#?# Again, what's this telling me.

#?# This looks like a further example?
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)




Please excuse the length of this email (and my lack of understanding)

Hope you can help and thanks.




Stephen Choularton Ph.D., FIoD


On 13/10/2010 2:41 AM, stefano iacus wrote:

just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:



Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn








------------------------------

Message: 2
Date: Tue, 12 Oct 2010 05:43:32 -0400
From: Sarbo
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
Message-ID:
Content-Type: text/plain; charset="utf-8"

By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to
discretise the equation and then fit the data to it. I don't have the
time right now to go into all the details but it's not hard- you can
parameterise the process using simple linear regression. If you need
help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:



Hi

Wonder if anyone could point me how I use this method to discover the
half life of a mean reverting process.

I am looking into pair trading and the time it takes for a
cointegrated pair to revert to the norm.

--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions


should go.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html


-------------- next part --------------
A non-text attachment was scrubbed...
Name: CoS2010Winner.JPG
Type: image/jpeg
Size: 16091 bytes
Desc: not available
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe


------------------------------

_______________________________________________
R-SIG-Finance mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-finance


End of R-SIG-Finance Digest, Vol 77, Issue 8
********************************************


[[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.



-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com







_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com




_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
Subscriber-posting only. If you want to post, subscribe first. -- Also
note that this is not the r-help list where general R questions should
go.
 		 	   		  
_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
No virus found in this incoming message. Checked by AVG - www.avg.com

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Mean reversion

Yihao Lu aeolus_lu

regarding your problem, i think Bjorn's link is very helpful. it even has the matlab code attached. several lines, you can translate it into R and do a comparison on the data, you will have the idea which theta is the rate you should look at.



Best,
Yihao







________________________________

> Date: Tue, 19 Oct 2010 12:50:46 +1100
> From: [hidden email]
> To: [hidden email]
> CC: [hidden email]
> Subject: Re: [R-SIG-Finance] Mean reversion
>
> Wish I could, but as you can see I am having even greater problems ;-)
>
> Stephen Choularton Ph.D., FIoD
>
> On 19/10/2010 12:35 PM, Yihao Lu aeolus_lu wrote:
>
>
> I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
> Is there anyone who can give me some possible explanation or guide me to some reference? thanks
>
> Best,
> Yihao
>
>
>
>
>
>
>
> ________________________________
>
>
> Date: Tue, 19 Oct 2010 09:03:55 +1100
> From: [hidden email]
> To: [hidden email]
> CC: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>
> Hi
>
> I am still trying to sort this one out. Any comments from anyone would
> be most welcome.
>
> Stephen Choularton Ph.D., FIoD
>
>
>
> On 14/10/2010 7:29 AM, Stephen Choularton wrote:
> Thanks for this help.
>
> Trying to make sense of it so I have added some notes to the code. I
> have marked them #?#
>
> Delighted if you can tell me if I am write or wrong, add any comments,
> answers.
>
>
> #?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
> #?# process work' particularly via dcOU
> #?# I have noted in several places that I am after:
> #?# 'the half-life of the decay equals ln(2)/θ'
> #?# 'The half-life is given as log(2)/mean-reversion speed.'
> #?# and I see theta appearing at a number of points in the code.
> #?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
> #?# eg is one of these the theta I am after?
>
> # ex3.01.R
> OU.lik <- function(theta1, theta2, theta3){
> n <- length(X)
> dt <- deltat(X)
> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> }
>
> require(stats4)
> require(sde)
>
> #?# random numer generation seed
> set.seed(123)
>
> #?# creation of a data set
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> #?# If I Look at X its like this:
> #?# Time Series:
> #?# Start = 0
> #?# End = 1000
> #?# Frequency = 1
> #?# [1] 1.00000000 etc
> #?# What sort of data object is it and how would I coerce an object with one
> #?# column from a read.csv into it?
>
>
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> summary(fit)
>
> #?# This gives:
>
> #?# Maximum likelihood estimation
>
> #?# Call:
> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>
> #?# Coefficients:
> #?# Estimate Std. Error
> #?# theta1 3.355322 0.28159504
> #?# theta2 1.106107 0.09010627
> #?# theta3 2.052815 0.07624441
>
> #?# -2 log L: 3366.389
>
> #?# What's this telling me?
>
> # ex3.01.R (cont.)
> prof <- profile(fit)
> par(mfrow=c(1,3))
> plot(prof)
> par(mfrow=c(1,1))
> vcov(fit)
> confint(fit)
>
> #?# This provides me with this output using 'fit' from before:
>
> #?# > vcov(fit)
> #?# theta1 theta2 theta3
> #?# theta1 0.07929576 0.024620718 0.016634557
> #?# theta2 0.02462072 0.008119141 0.005485549
> #?# theta3 0.01663456 0.005485549 0.005813209
> #?# > confint(fit)
> #?# Profiling...
> #?# 2.5 % 97.5 %
> #?# theta1 2.8448980 3.960982
> #?# theta2 0.9433338 1.300629
> #?# theta3 1.9147136 2.216113
>
> #?# and 'fit' is:
>
> #?# Call:
> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>
> #?# Coefficients:
> #?# theta1 theta2 theta3
> #?# 3.355322 1.106107 2.052815
>
> #?# plus some graphic output
>
> #?# Again, what's this telling me.
>
> #?# This looks like a further example?
> # ex3.01.R (cont.)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> summary(fit2)
>
>
>
>
> Please excuse the length of this email (and my lack of understanding)
>
> Hope you can help and thanks.
>
>
>
>
> Stephen Choularton Ph.D., FIoD
>
>
> On 13/10/2010 2:41 AM, stefano iacus wrote:
>
> just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
>
> sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
>
> This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
>
>
> # ex3.01.R
> OU.lik <- function(theta1, theta2, theta3){
> n <- length(X)
> dt <- deltat(X)
> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> }
>
> require(stats4)
> require(sde)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> summary(fit)
>
> # ex3.01.R (cont.)
> prof <- profile(fit)
> par(mfrow=c(1,3))
> plot(prof)
> par(mfrow=c(1,1))
> vcov(fit)
> confint(fit)
>
> # ex3.01.R (cont.)
> set.seed(123)
> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> summary(fit2)
>
>
> I hope this helps out
>
> stefano
>
> On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
>
>
>
> Hi Stephen,
>
> You could take a look at
>
> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>
> for the linear regression method, or take a look at the package "sde" which
> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
> though, only the CIR).
>
> The half-life is given as log(2)/mean-reversion speed.
>
> Do keep an eye on the partition of the time-axis, e.g. what frequency you
> are using (daily, yearly) for interpreting the half-life.
>
> BR,
> Bjørn
>
>
>
>
>
>
>
>
> ------------------------------
>
> Message: 2
> Date: Tue, 12 Oct 2010 05:43:32 -0400
> From: Sarbo
> To: [hidden email]
> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> Message-ID:
> Content-Type: text/plain; charset="utf-8"
>
> By half-life, do you mean the speed of mean-reversion?
>
> If so, there's a bit of algebraic tomfoolery that's required to
> discretise the equation and then fit the data to it. I don't have the
> time right now to go into all the details but it's not hard- you can
> parameterise the process using simple linear regression. If you need
> help with that I'll try and get back to you tonight about it.
>
> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>
>
>
> Hi
>
> Wonder if anyone could point me how I use this method to discover the
> half life of a mean reverting process.
>
> I am looking into pair trading and the time it takes for a
> cointegrated pair to revert to the norm.
>
> --
> Stephen Choularton Ph.D., FIoD
>
> 9999 2226
> 0413 545 182
>
>
> for insurance go to www.netinsure.com.au
> for markets go to www.organicfoodmarkets.com.au
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions
>
>
> should go.
>
>
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>
>
> -------------- next part --------------
> A non-text attachment was scrubbed...
> Name: CoS2010Winner.JPG
> Type: image/jpeg
> Size: 16091 bytes
> Desc: not available
> URL: <
> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>
>
> ------------------------------
>
> _______________________________________________
> R-SIG-Finance mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>
>
> End of R-SIG-Finance Digest, Vol 77, Issue 8
> ********************************************
>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
> -----------------------------------
> Stefano M. Iacus
> Department of Economics,
> Business and Statistics
> University of Milan
> Via Conservatorio, 7
> I-20123 Milan - Italy
> Ph.: +39 02 50321 461
> Fax: +39 02 50321 505
> http://www.economia.unimi.it/iacus
> ------------------------------------------------------------------------------------
> Please don't send me Word or PowerPoint attachments if not
> absolutely necessary. See:
> http://www.gnu.org/philosophy/no-word-attachments.html
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
>
>
>
>
>
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
>
>
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
> Subscriber-posting only. If you want to post, subscribe first. -- Also
> note that this is not the r-help list where general R questions should
> go.
>
>
>
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.
>
>
>
>
>
> No virus found in this incoming message.
> Checked by AVG - www.avg.com
> Version: 9.0.862 / Virus Database: 271.1.1/3204 - Release Date: 10/18/10 17:34:00
>
>

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

cointegration

Stephen Choularton-3
In reply to this post by Yihao Lu aeolus_lu
Hi Folks

I'm using this to find cointegrated stocks on the AX.

library(xts)
library(quantmod)

# quickly re-source this file
s <- function() source('meanrev.R')

checkPairFromYahoo <- function(sym1, sym2, dateFilter='::')
{
  t.xts <- getCombined(sym1, sym2, dateFilter=dateFilter)

  cat("Date range is", format(start(t.xts)), "to", format(end(t.xts)), "\n")

  # Build linear model
  m <- buildLM(t.xts)

  # Note beta -- http://en.wikipedia.org/wiki/Beta_(finance)
  beta <- getBeta(m)
  cat("Assumed hedge ratio is", beta, "\n")

  # Build spread
  sprd <- buildSpread(t.xts, beta)

  # Test cointegration
  ht <- testCoint(sprd)
  cat("PP p-value is", as.double(ht$p.value), "\n")

  if (as.double(ht$p.value) < 0.05)
  {
    cat("###############################################################\n", sym1 ,":", sym2 ," is likely mean-reverting.\n", "###########################################################\n" )
  }
  else
  {
    #cat(sym1 ,":", sym2 ," is not mean-reverting.\n")
  }
}

getCombined <- function(sym1, sym2, dateFilter='::')
{
  # Grab historical data for both symbols
  one <- getSymbols(sym1, auto.assign=FALSE)
  two <- getSymbols(sym2, auto.assign=FALSE)

  # Give columns more usable names
  colnames(one) <- c('Open', 'High', 'Low', 'Close', 'Volume', 'Adjusted')
  colnames(two) <- c('Open', 'High', 'Low', 'Close', 'Volume', 'Adjusted')

  # Build combined object
  return(merge(one$Close, two$Close, all=FALSE)[dateFilter])
}

buildLM <- function(combined)
{
  return(lm(Close ~ Close.1 + 0, combined))
}

getBeta <- function(m)
{
  return(as.double(coef(m)[1]))
}

buildSpread <- function(combined, beta)
{
  return(combined$Close - beta*combined$Close.1)
}

testCoint <- function(sprd)
{
  return(PP.test(sprd, lshort = FALSE))
}

I run it on batches of stock-pairs and then have a look at those which are cointegrated.  Assuming my code is right (and anyone who thinks there is something wrong with it please let me know ;-)

Just wondered if anyone simply goes with the results, or if a test of logic is required.  I found, for example, that AGL ( a big gas company) was cointegrated with Bunnings Wharehouses (a hardware superstore chain).  Can't see the reason for that.  AMP (major insurer) cointegrates with AXA (another major insurer).  That makes sense and it cointegrates with  Westpac (major bank) still some logic but a bit thinner.  It also cointegrates with Fortescue Metals (big iron ore operation).  Not much logic there.  Anyway question is: do you get better results by using informed judgement on these things or just trust the figures?

Any comments most welcome.


Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au

On 19/10/2010 12:35 PM, Yihao Lu aeolus_lu wrote:
I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
Is there anyone who can give me some possible explanation or guide me to some reference? thanks

Best,
Yihao







________________________________
Date: Tue, 19 Oct 2010 09:03:55 +1100
From: [hidden email]
To: [hidden email]
CC: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck

Hi

I am still trying to sort this one out. Any comments from anyone would
be most welcome.

Stephen Choularton Ph.D., FIoD



On 14/10/2010 7:29 AM, Stephen Choularton wrote:
Thanks for this help.

Trying to make sense of it so I have added some notes to the code. I
have marked them #?#

Delighted if you can tell me if I am write or wrong, add any comments,
answers.


#?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
#?# process work' particularly via dcOU
#?# I have noted in several places that I am after:
#?# 'the half-life of the decay equals ln(2)/θ'
#?# 'The half-life is given as log(2)/mean-reversion speed.'
#?# and I see theta appearing at a number of points in the code.
#?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
#?# eg is one of these the theta I am after?

# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)

#?# random numer generation seed
set.seed(123)

#?# creation of a data set
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
#?# If I Look at X its like this:
#?# Time Series:
#?# Start = 0
#?# End = 1000
#?# Frequency = 1
#?# [1] 1.00000000 etc
#?# What sort of data object is it and how would I coerce an object with one
#?# column from a read.csv into it?


mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

#?# This gives:

#?# Maximum likelihood estimation

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# Estimate Std. Error
#?# theta1 3.355322 0.28159504
#?# theta2 1.106107 0.09010627
#?# theta3 2.052815 0.07624441

#?# -2 log L: 3366.389

#?# What's this telling me?

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

#?# This provides me with this output using 'fit' from before:

#?# > vcov(fit)
#?# theta1 theta2 theta3
#?# theta1 0.07929576 0.024620718 0.016634557
#?# theta2 0.02462072 0.008119141 0.005485549
#?# theta3 0.01663456 0.005485549 0.005813209
#?# > confint(fit)
#?# Profiling...
#?# 2.5 % 97.5 %
#?# theta1 2.8448980 3.960982
#?# theta2 0.9433338 1.300629
#?# theta3 1.9147136 2.216113

#?# and 'fit' is:

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# theta1 theta2 theta3
#?# 3.355322 1.106107 2.052815

#?# plus some graphic output

#?# Again, what's this telling me.

#?# This looks like a further example?
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)




Please excuse the length of this email (and my lack of understanding)

Hope you can help and thanks.




Stephen Choularton Ph.D., FIoD


On 13/10/2010 2:41 AM, stefano iacus wrote:

just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:



Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn








------------------------------

Message: 2
Date: Tue, 12 Oct 2010 05:43:32 -0400
From: Sarbo
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
Message-ID:
Content-Type: text/plain; charset="utf-8"

By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to
discretise the equation and then fit the data to it. I don't have the
time right now to go into all the details but it's not hard- you can
parameterise the process using simple linear regression. If you need
help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:



Hi

Wonder if anyone could point me how I use this method to discover the
half life of a mean reverting process.

I am looking into pair trading and the time it takes for a
cointegrated pair to revert to the norm.

--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions


should go.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html


-------------- next part --------------
A non-text attachment was scrubbed...
Name: CoS2010Winner.JPG
Type: image/jpeg
Size: 16091 bytes
Desc: not available
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe


------------------------------

_______________________________________________
R-SIG-Finance mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-finance


End of R-SIG-Finance Digest, Vol 77, Issue 8
********************************************


[[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.



-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com







_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com




_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
Subscriber-posting only. If you want to post, subscribe first. -- Also
note that this is not the r-help list where general R questions should
go.
 		 	   		  
_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
No virus found in this incoming message. Checked by AVG - www.avg.com

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Mean reversion

Eric Zivot
In reply to this post by Yihao Lu aeolus_lu
You have to be careful when interpreting rolling ADF tests because the usual critical values for evaluating the tests are not valid. See the 1992 Journal of Business and Economic Statistics paper by Banerjee, Lumsdaine and Stock.

****************************************************************
*  Eric Zivot                                 *
*  Professor and Gary Waterman Distinguished Scholar           *
*  Department of Economics                                     *
*  Adjunct Professor of Finance                                *
*  Adjunct Professor of Statistics
*  Box 353330                  email:  [hidden email] *
*  University of Washington    phone:  206-543-6715            *
*  Seattle, WA 98195-3330                                      *                                                           *
*  www:  http://faculty.washington.edu/ezivot                  *
****************************************************************

On Mon, 18 Oct 2010, Yihao Lu aeolus_lu wrote:

>
> I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
> Is there anyone who can give me some possible explanation or guide me to some reference? thanks
>
> Best,
> Yihao
>
>
>
>
>
>
>
> ________________________________
>> Date: Tue, 19 Oct 2010 09:03:55 +1100
>> From: [hidden email]
>> To: [hidden email]
>> CC: [hidden email]
>> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>>
>> Hi
>>
>> I am still trying to sort this one out. Any comments from anyone would
>> be most welcome.
>>
>> Stephen Choularton Ph.D., FIoD
>>
>>
>>
>> On 14/10/2010 7:29 AM, Stephen Choularton wrote:
>> Thanks for this help.
>>
>> Trying to make sense of it so I have added some notes to the code. I
>> have marked them #?#
>>
>> Delighted if you can tell me if I am write or wrong, add any comments,
>> answers.
>>
>>
>> #?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
>> #?# process work' particularly via dcOU
>> #?# I have noted in several places that I am after:
>> #?# 'the half-life of the decay equals ln(2)/θ'
>> #?# 'The half-life is given as log(2)/mean-reversion speed.'
>> #?# and I see theta appearing at a number of points in the code.
>> #?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
>> #?# eg is one of these the theta I am after?
>>
>> # ex3.01.R
>> OU.lik <- function(theta1, theta2, theta3){
>> n <- length(X)
>> dt <- deltat(X)
>> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
>> }
>>
>> require(stats4)
>> require(sde)
>>
>> #?# random numer generation seed
>> set.seed(123)
>>
>> #?# creation of a data set
>> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
>> #?# If I Look at X its like this:
>> #?# Time Series:
>> #?# Start = 0
>> #?# End = 1000
>> #?# Frequency = 1
>> #?# [1] 1.00000000 etc
>> #?# What sort of data object is it and how would I coerce an object with one
>> #?# column from a read.csv into it?
>>
>>
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
>> summary(fit)
>>
>> #?# This gives:
>>
>> #?# Maximum likelihood estimation
>>
>> #?# Call:
>> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
>> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>>
>> #?# Coefficients:
>> #?# Estimate Std. Error
>> #?# theta1 3.355322 0.28159504
>> #?# theta2 1.106107 0.09010627
>> #?# theta3 2.052815 0.07624441
>>
>> #?# -2 log L: 3366.389
>>
>> #?# What's this telling me?
>>
>> # ex3.01.R (cont.)
>> prof <- profile(fit)
>> par(mfrow=c(1,3))
>> plot(prof)
>> par(mfrow=c(1,1))
>> vcov(fit)
>> confint(fit)
>>
>> #?# This provides me with this output using 'fit' from before:
>>
>> #?# > vcov(fit)
>> #?# theta1 theta2 theta3
>> #?# theta1 0.07929576 0.024620718 0.016634557
>> #?# theta2 0.02462072 0.008119141 0.005485549
>> #?# theta3 0.01663456 0.005485549 0.005813209
>> #?# > confint(fit)
>> #?# Profiling...
>> #?# 2.5 % 97.5 %
>> #?# theta1 2.8448980 3.960982
>> #?# theta2 0.9433338 1.300629
>> #?# theta3 1.9147136 2.216113
>>
>> #?# and 'fit' is:
>>
>> #?# Call:
>> #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
>> #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
>>
>> #?# Coefficients:
>> #?# theta1 theta2 theta3
>> #?# 3.355322 1.106107 2.052815
>>
>> #?# plus some graphic output
>>
>> #?# Again, what's this telling me.
>>
>> #?# This looks like a further example?
>> # ex3.01.R (cont.)
>> set.seed(123)
>> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
>> summary(fit2)
>>
>>
>>
>>
>> Please excuse the length of this email (and my lack of understanding)
>>
>> Hope you can help and thanks.
>>
>>
>>
>>
>> Stephen Choularton Ph.D., FIoD
>>
>>
>> On 13/10/2010 2:41 AM, stefano iacus wrote:
>>
>> just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
>>
>> sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
>>
>> This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
>>
>>
>> # ex3.01.R
>> OU.lik <- function(theta1, theta2, theta3){
>> n <- length(X)
>> dt <- deltat(X)
>> -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
>> }
>>
>> require(stats4)
>> require(sde)
>> set.seed(123)
>> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
>> summary(fit)
>>
>> # ex3.01.R (cont.)
>> prof <- profile(fit)
>> par(mfrow=c(1,3))
>> plot(prof)
>> par(mfrow=c(1,1))
>> vcov(fit)
>> confint(fit)
>>
>> # ex3.01.R (cont.)
>> set.seed(123)
>> X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
>> mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
>> method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
>> summary(fit2)
>>
>>
>> I hope this helps out
>>
>> stefano
>>
>> On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
>>
>>
>>
>> Hi Stephen,
>>
>> You could take a look at
>>
>> http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
>>
>> for the linear regression method, or take a look at the package "sde" which
>> contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
>> though, only the CIR).
>>
>> The half-life is given as log(2)/mean-reversion speed.
>>
>> Do keep an eye on the partition of the time-axis, e.g. what frequency you
>> are using (daily, yearly) for interpreting the half-life.
>>
>> BR,
>> Bjørn
>>
>>
>>
>>
>>
>>
>>
>>
>> ------------------------------
>>
>> Message: 2
>> Date: Tue, 12 Oct 2010 05:43:32 -0400
>> From: Sarbo
>> To: [hidden email]
>> Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
>> Message-ID:
>> Content-Type: text/plain; charset="utf-8"
>>
>> By half-life, do you mean the speed of mean-reversion?
>>
>> If so, there's a bit of algebraic tomfoolery that's required to
>> discretise the equation and then fit the data to it. I don't have the
>> time right now to go into all the details but it's not hard- you can
>> parameterise the process using simple linear regression. If you need
>> help with that I'll try and get back to you tonight about it.
>>
>> On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
>>
>>
>>
>> Hi
>>
>> Wonder if anyone could point me how I use this method to discover the
>> half life of a mean reverting process.
>>
>> I am looking into pair trading and the time it takes for a
>> cointegrated pair to revert to the norm.
>>
>> --
>> Stephen Choularton Ph.D., FIoD
>>
>> 9999 2226
>> 0413 545 182
>>
>>
>> for insurance go to www.netinsure.com.au
>> for markets go to www.organicfoodmarkets.com.au
>>
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions
>>
>>
>> should go.
>>
>>
>> -------------- next part --------------
>> An HTML attachment was scrubbed...
>> URL: <
>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
>>
>>
>> -------------- next part --------------
>> A non-text attachment was scrubbed...
>> Name: CoS2010Winner.JPG
>> Type: image/jpeg
>> Size: 16091 bytes
>> Desc: not available
>> URL: <
>> https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
>>
>>
>> ------------------------------
>>
>> _______________________________________________
>> R-SIG-Finance mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>>
>>
>> End of R-SIG-Finance Digest, Vol 77, Issue 8
>> ********************************************
>>
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions should go.
>>
>>
>>
>> -----------------------------------
>> Stefano M. Iacus
>> Department of Economics,
>> Business and Statistics
>> University of Milan
>> Via Conservatorio, 7
>> I-20123 Milan - Italy
>> Ph.: +39 02 50321 461
>> Fax: +39 02 50321 505
>> http://www.economia.unimi.it/iacus
>> ------------------------------------------------------------------------------------
>> Please don't send me Word or PowerPoint attachments if not
>> absolutely necessary. See:
>> http://www.gnu.org/philosophy/no-word-attachments.html
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions should go.
>>
>>
>>
>>
>>
>> No virus found in this incoming message.
>> Checked by AVG - www.avg.com
>>
>>
>>
>>
>>
>>
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
>> -- Subscriber-posting only. If you want to post, subscribe first.
>> -- Also note that this is not the r-help list where general R questions should go.
>>
>>
>>
>>
>>
>> No virus found in this incoming message.
>> Checked by AVG - www.avg.com
>>
>>
>>
>>
>> _______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
>> Subscriber-posting only. If you want to post, subscribe first. -- Also
>> note that this is not the r-help list where general R questions should
>> go.
>    
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Mean reversion

Sarbo
In reply to this post by Yihao Lu aeolus_lu
Yihao,

Prof. Zivot is right. The ADF test isn't a great way to test for
mean-reversion; it's merely a way to test for stationarity to specify
the degrees in an ARIMA model.

I dug up some code from my misspent youth as a consultant which you
might find useful:

StochasticProcessTest <- function(PriceSeries, delta.t, p.value,
diagnostics = TRUE, from = NULL, to = NULL,
                                  by = NULL, currency = '($)'){
  require(stats)
  require(fSeries)

  m <- length(PriceSeries)
  S1 <- PriceSeries[1:(m-1)]
  S2 <- PriceSeries[2:m]

  #Fit a GBM process:
  Y <- diff(log(PriceSeries)); X <- S1
  GBMfit <- lm(Y ~ X)

  #Fit a mean-reverting GBM process:
  Y <- diff(PriceSeries) / S1; X <- log(S1)
  MRGBMfit <- lm(Y ~ X)

  #Fit a Vasicek (O-U) process:
  Y <- diff(PriceSeries); X <- S1
  Vfit <- lm(Y ~ X)

  #Fit a Cox-Ingersoll-Ross process:
  Y <- diff(PriceSeries) / sqrt(S1); X <- S1 / sqrt(S1)
  CIRfit <- lm(Y ~ X)

  #Gather the parameter estimates:
  kappa <- as.vector((c(0, -MRGBMfit$coef[2], -Vfit$coef[2], -CIRfit
$coef[2])) / delta.t)
  mu <- as.vector((c(mean(log(S2/S1)), MRGBMfit$coef[1] / kappa[2], Vfit
$coef[1] / kappa[3],
                     CIRfit$coef[1] / kappa[4])) / delta.t)
  sigma <- (c(sd(GBMfit$resid), sd(MRGBMfit$resid), sd(Vfit$resid),
sd(CIRfit$resid))) / sqrt(delta.t)
  tstat <- as.vector(c(summary(GBMfit)$coef[2,3],
summary(MRGBMfit)$coef[2,3], summary(Vfit)$coef[2,3],
                       summary(CIRfit)$coef[2,3]))
  fstat <- list(summary(GBMfit)$fstatistic,
summary(MRGBMfit)$fstatistic, summary(Vfit)$fstatistic,
                summary(CIRfit)$fstatistic)
  names(fstat) <- c('GBM', 'Mean.Revert.GBM', 'Vasicek', 'CIR')
  AICs <- c(AIC(GBMfit), AIC(MRGBMfit), AIC(Vfit), AIC(CIRfit))
  paramframe <- data.frame(rbind(kappa, mu, sigma, tstat, AICs))
  names(paramframe) <- names(fstat)

  #Now figure out what the actual process is, using AIC:
  crit <- ifelse(m > 30, qnorm(1 - p.value), qt(1 - p.value, df = n))
  tmp <- which.min(AICs)
  Processes <- names(fstat)
  Verdict <- Processes[tmp]
  FinalSummary <- switch(Verdict, GBM = summary(GBMfit), CIR =
summary(CIRfit), Mean.Revert.GBM = MRGBMfit,
                         Vasicek = Vfit)
  fitobj <- switch(Verdict, GBM = GBMfit, CIR = CIRfit, Mean.Revert.GBM
= MRGBMfit, Vasicek = Vfit)
  Output <- list(Parameters = paramframe, Critical.Value = crit, Verdict
= Verdict, FinalSummary = FinalSummary,
                 fstat = fstat, fitted.object = fitobj)

  if (diagnostics){
    op <- par(ask = TRUE)
    on.exit(op)

    if(!all(c(class(from), class(to)) == 'Date')){
      S <- timeSeries(PriceSeries)
    } else S = timeSeries(PriceSeries, seq(from, to, length.out = m))

    plot(S, type = 'l', xlab = 'Date', ylab = paste('Price Series',
currency), main = 'Time Series Plot of Data',
         lwd = 2, col = 'blue')
    rets <- returns(PriceSeries, 'continuous')[-1]
    hist(rets, xlab = paste('Log-Returns', currency), col = 'blue',
border = 'white', main = 'Histogram of Return Series',
         freq = FALSE)
    x <- seq(min(rets), max(rets), length.out = max(m, 1000))
    lines(x, dnorm(x, mean(rets), sd(rets)), col = 'magenta', lwd = 2)
    lines(density(rets), col = 'green', lwd = 2)
    legend('topright', legend = c('Log-Returns', 'Observed CDF',
'Gaussian Fit'), lwd = rep(2, 3),
                         col = c('blue', 'green', 'magenta'))
    plot(fitobj, 1:6)
  }

  return(Output)

}
   

(I don't claim that it's necessarily great code, but it does seem to
work.)

On Mon, 2010-10-18 at 21:35 -0400, Yihao Lu aeolus_lu wrote:

> I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
> Is there anyone who can give me some possible explanation or guide me to some reference? thanks
>
> Best,
> Yihao
>
>
>
>
>
>
>
> ________________________________
> > Date: Tue, 19 Oct 2010 09:03:55 +1100
> > From: [hidden email]
> > To: [hidden email]
> > CC: [hidden email]
> > Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> >
> > Hi
> >
> > I am still trying to sort this one out. Any comments from anyone would
> > be most welcome.
> >
> > Stephen Choularton Ph.D., FIoD
> >
> >
> >
> > On 14/10/2010 7:29 AM, Stephen Choularton wrote:
> > Thanks for this help.
> >
> > Trying to make sense of it so I have added some notes to the code. I
> > have marked them #?#
> >
> > Delighted if you can tell me if I am write or wrong, add any comments,
> > answers.
> >
> >
> > #?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
> > #?# process work' particularly via dcOU
> > #?# I have noted in several places that I am after:
> > #?# 'the half-life of the decay equals ln(2)/θ'
> > #?# 'The half-life is given as log(2)/mean-reversion speed.'
> > #?# and I see theta appearing at a number of points in the code.
> > #?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
> > #?# eg is one of these the theta I am after?
> >
> > # ex3.01.R
> > OU.lik <- function(theta1, theta2, theta3){
> > n <- length(X)
> > dt <- deltat(X)
> > -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> > }
> >
> > require(stats4)
> > require(sde)
> >
> > #?# random numer generation seed
> > set.seed(123)
> >
> > #?# creation of a data set
> > X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> > #?# If I Look at X its like this:
> > #?# Time Series:
> > #?# Start = 0
> > #?# End = 1000
> > #?# Frequency = 1
> > #?# [1] 1.00000000 etc
> > #?# What sort of data object is it and how would I coerce an object with one
> > #?# column from a read.csv into it?
> >
> >
> > mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> > method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> > summary(fit)
> >
> > #?# This gives:
> >
> > #?# Maximum likelihood estimation
> >
> > #?# Call:
> > #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> > #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
> >
> > #?# Coefficients:
> > #?# Estimate Std. Error
> > #?# theta1 3.355322 0.28159504
> > #?# theta2 1.106107 0.09010627
> > #?# theta3 2.052815 0.07624441
> >
> > #?# -2 log L: 3366.389
> >
> > #?# What's this telling me?
> >
> > # ex3.01.R (cont.)
> > prof <- profile(fit)
> > par(mfrow=c(1,3))
> > plot(prof)
> > par(mfrow=c(1,1))
> > vcov(fit)
> > confint(fit)
> >
> > #?# This provides me with this output using 'fit' from before:
> >
> > #?# > vcov(fit)
> > #?# theta1 theta2 theta3
> > #?# theta1 0.07929576 0.024620718 0.016634557
> > #?# theta2 0.02462072 0.008119141 0.005485549
> > #?# theta3 0.01663456 0.005485549 0.005813209
> > #?# > confint(fit)
> > #?# Profiling...
> > #?# 2.5 % 97.5 %
> > #?# theta1 2.8448980 3.960982
> > #?# theta2 0.9433338 1.300629
> > #?# theta3 1.9147136 2.216113
> >
> > #?# and 'fit' is:
> >
> > #?# Call:
> > #?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
> > #?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))
> >
> > #?# Coefficients:
> > #?# theta1 theta2 theta3
> > #?# 3.355322 1.106107 2.052815
> >
> > #?# plus some graphic output
> >
> > #?# Again, what's this telling me.
> >
> > #?# This looks like a further example?
> > # ex3.01.R (cont.)
> > set.seed(123)
> > X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> > mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> > method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> > summary(fit2)
> >
> >
> >
> >
> > Please excuse the length of this email (and my lack of understanding)
> >
> > Hope you can help and thanks.
> >
> >
> >
> >
> > Stephen Choularton Ph.D., FIoD
> >
> >
> > On 13/10/2010 2:41 AM, stefano iacus wrote:
> >
> > just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.
> >
> > sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.
> >
> > This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU
> >
> >
> > # ex3.01.R
> > OU.lik <- function(theta1, theta2, theta3){
> > n <- length(X)
> > dt <- deltat(X)
> > -sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
> > }
> >
> > require(stats4)
> > require(sde)
> > set.seed(123)
> > X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
> > mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> > method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
> > summary(fit)
> >
> > # ex3.01.R (cont.)
> > prof <- profile(fit)
> > par(mfrow=c(1,3))
> > plot(prof)
> > par(mfrow=c(1,1))
> > vcov(fit)
> > confint(fit)
> >
> > # ex3.01.R (cont.)
> > set.seed(123)
> > X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
> > mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
> > method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
> > summary(fit2)
> >
> >
> > I hope this helps out
> >
> > stefano
> >
> > On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:
> >
> >
> >
> > Hi Stephen,
> >
> > You could take a look at
> >
> > http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model
> >
> > for the linear regression method, or take a look at the package "sde" which
> > contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
> > though, only the CIR).
> >
> > The half-life is given as log(2)/mean-reversion speed.
> >
> > Do keep an eye on the partition of the time-axis, e.g. what frequency you
> > are using (daily, yearly) for interpreting the half-life.
> >
> > BR,
> > Bjørn
> >
> >
> >
> >
> >
> >
> >
> >
> > ------------------------------
> >
> > Message: 2
> > Date: Tue, 12 Oct 2010 05:43:32 -0400
> > From: Sarbo
> > To: [hidden email]
> > Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
> > Message-ID:
> > Content-Type: text/plain; charset="utf-8"
> >
> > By half-life, do you mean the speed of mean-reversion?
> >
> > If so, there's a bit of algebraic tomfoolery that's required to
> > discretise the equation and then fit the data to it. I don't have the
> > time right now to go into all the details but it's not hard- you can
> > parameterise the process using simple linear regression. If you need
> > help with that I'll try and get back to you tonight about it.
> >
> > On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:
> >
> >
> >
> > Hi
> >
> > Wonder if anyone could point me how I use this method to discover the
> > half life of a mean reverting process.
> >
> > I am looking into pair trading and the time it takes for a
> > cointegrated pair to revert to the norm.
> >
> > --
> > Stephen Choularton Ph.D., FIoD
> >
> > 9999 2226
> > 0413 545 182
> >
> >
> > for insurance go to www.netinsure.com.au
> > for markets go to www.organicfoodmarkets.com.au
> >
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R questions
> >
> >
> > should go.
> >
> >
> > -------------- next part --------------
> > An HTML attachment was scrubbed...
> > URL: <
> > https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html
> >
> >
> > -------------- next part --------------
> > A non-text attachment was scrubbed...
> > Name: CoS2010Winner.JPG
> > Type: image/jpeg
> > Size: 16091 bytes
> > Desc: not available
> > URL: <
> > https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe
> >
> >
> > ------------------------------
> >
> > _______________________________________________
> > R-SIG-Finance mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> >
> >
> > End of R-SIG-Finance Digest, Vol 77, Issue 8
> > ********************************************
> >
> >
> > [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R questions should go.
> >
> >
> >
> > -----------------------------------
> > Stefano M. Iacus
> > Department of Economics,
> > Business and Statistics
> > University of Milan
> > Via Conservatorio, 7
> > I-20123 Milan - Italy
> > Ph.: +39 02 50321 461
> > Fax: +39 02 50321 505
> > http://www.economia.unimi.it/iacus
> > ------------------------------------------------------------------------------------
> > Please don't send me Word or PowerPoint attachments if not
> > absolutely necessary. See:
> > http://www.gnu.org/philosophy/no-word-attachments.html
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R questions should go.
> >
> >
> >
> >
> >
> > No virus found in this incoming message.
> > Checked by AVG - www.avg.com
> >
> >
> >
> >
> >
> >
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> > -- Subscriber-posting only. If you want to post, subscribe first.
> > -- Also note that this is not the r-help list where general R questions should go.
> >
> >
> >
> >
> >
> > No virus found in this incoming message.
> > Checked by AVG - www.avg.com
> >
> >
> >
> >
> > _______________________________________________
> > [hidden email] mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
> > Subscriber-posting only. If you want to post, subscribe first. -- Also
> > note that this is not the r-help list where general R questions should
> > go.
>      
> _______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-finance
> -- Subscriber-posting only. If you want to post, subscribe first.
> -- Also note that this is not the r-help list where general R questions should go.


        [[alternative HTML version deleted]]


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cointegration

Paul Teetor
In reply to this post by Stephen Choularton-3

Stephen,
 
It depends what you mean by "logic".
 
If you mean statistical logic, I'll defer to Eric Zivot and Sarbo who are far wiser than I am. I will note, however, that you are testing for a p-value of 0.05, so I expect 5% of your test results to be misleading. In other words, for every 20 pairs tested by your batch job, I expect one will be suspect.
 
"Spurious cointegration" is a serious problem. I suggest Googling that topic. You may be suprised what you learn. (The irony, of course, is that cointegration was supposed to cure "spurious correlation." Oh well.)
 
If you mean financial logic, I strongly suggest not blindly risking money on your statistical test. Some filtering is required. Look for trades that make sense.
 
For example, my software reports that the stocks of MSFT and GOOG form a mean-reverting pair. But I would not trade that spread: too much idiosyncratic risk. My software also reports that Corn futures and Soybean Oil futures form a mean-reverting pair. But I would not trade that spread because the economic connection between corn and bean oil is too weak.
 
Hope that helps.
 
Paul


From: [hidden email] [mailto:[hidden email]] On Behalf Of Stephen Choularton
Sent: Monday, October 18, 2010 9:46 PM
To: [hidden email]
Subject: [R-SIG-Finance] cointegration

Hi Folks

I'm using this to find cointegrated stocks on the AX.

library(xts)
library(quantmod)

# quickly re-source this file
s <- function() source('meanrev.R')

checkPairFromYahoo <- function(sym1, sym2, dateFilter='::')
{
  t.xts <- getCombined(sym1, sym2, dateFilter=dateFilter)

  cat("Date range is", format(start(t.xts)), "to", format(end(t.xts)), "\n")

  # Build linear model
  m <- buildLM(t.xts)

  # Note beta -- http://en.wikipedia.org/wiki/Beta_(finance)
  beta <- getBeta(m)
  cat("Assumed hedge ratio is", beta, "\n")

  # Build spread
  sprd <- buildSpread(t.xts, beta)

  # Test cointegration
  ht <- testCoint(sprd)
  cat("PP p-value is", as.double(ht$p.value), "\n")

  if (as.double(ht$p.value) < 0.05)
  {
    cat("###############################################################\n", sym1 ,":", sym2 ," is likely mean-reverting.\n", "###########################################################\n" )
  }
  else
  {
    #cat(sym1 ,":", sym2 ," is not mean-reverting.\n")
  }
}

getCombined <- function(sym1, sym2, dateFilter='::')
{
  # Grab historical data for both symbols
  one <- getSymbols(sym1, auto.assign=FALSE)
  two <- getSymbols(sym2, auto.assign=FALSE)

  # Give columns more usable names
  colnames(one) <- c('Open', 'High', 'Low', 'Close', 'Volume', 'Adjusted')
  colnames(two) <- c('Open', 'High', 'Low', 'Close', 'Volume', 'Adjusted')

  # Build combined object
  return(merge(one$Close, two$Close, all=FALSE)[dateFilter])
}

buildLM <- function(combined)
{
  return(lm(Close ~ Close.1 + 0, combined))
}

getBeta <- function(m)
{
  return(as.double(coef(m)[1]))
}

buildSpread <- function(combined, beta)
{
  return(combined$Close - beta*combined$Close.1)
}

testCoint <- function(sprd)
{
  return(PP.test(sprd, lshort = FALSE))
}

I run it on batches of stock-pairs and then have a look at those which are cointegrated.  Assuming my code is right (and anyone who thinks there is something wrong with it please let me know ;-)

Just wondered if anyone simply goes with the results, or if a test of logic is required.  I found, for example, that AGL ( a big gas company) was cointegrated with Bunnings Wharehouses (a hardware superstore chain).  Can't see the reason for that.  AMP (major insurer) cointegrates with AXA (another major insurer).  That makes sense and it cointegrates with  Westpac (major bank) still some logic but a bit thinner.  It also cointegrates with Fortescue Metals (big iron ore operation).  Not much logic there.  Anyway question is: do you get better results by using informed judgement on these things or just trust the figures?

Any comments most welcome.


Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au

On 19/10/2010 12:35 PM, Yihao Lu aeolus_lu wrote:
I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
Is there anyone who can give me some possible explanation or guide me to some reference? thanks

Best,
Yihao







________________________________
Date: Tue, 19 Oct 2010 09:03:55 +1100
From: [hidden email]
To: [hidden email]
CC: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck

Hi

I am still trying to sort this one out. Any comments from anyone would
be most welcome.

Stephen Choularton Ph.D., FIoD



On 14/10/2010 7:29 AM, Stephen Choularton wrote:
Thanks for this help.

Trying to make sense of it so I have added some notes to the code. I
have marked them #?#

Delighted if you can tell me if I am write or wrong, add any comments,
answers.


#?# This appears to be the function that is doing the 'Ornstein-Uhlenbeck
#?# process work' particularly via dcOU
#?# I have noted in several places that I am after:
#?# 'the half-life of the decay equals ln(2)/θ'
#?# 'The half-life is given as log(2)/mean-reversion speed.'
#?# and I see theta appearing at a number of points in the code.
#?# Can you tell me why 3 thetas viz theta1, theta2, theta3 and what they do?
#?# eg is one of these the theta I am after?

# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)

#?# random numer generation seed
set.seed(123)

#?# creation of a data set
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
#?# If I Look at X its like this:
#?# Time Series:
#?# Start = 0
#?# End = 1000
#?# Frequency = 1
#?# [1] 1.00000000 etc
#?# What sort of data object is it and how would I coerce an object with one
#?# column from a read.csv into it?


mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

#?# This gives:

#?# Maximum likelihood estimation

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# Estimate Std. Error
#?# theta1 3.355322 0.28159504
#?# theta2 1.106107 0.09010627
#?# theta3 2.052815 0.07624441

#?# -2 log L: 3366.389

#?# What's this telling me?

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

#?# This provides me with this output using 'fit' from before:

#?# > vcov(fit)
#?# theta1 theta2 theta3
#?# theta1 0.07929576 0.024620718 0.016634557
#?# theta2 0.02462072 0.008119141 0.005485549
#?# theta3 0.01663456 0.005485549 0.005813209
#?# > confint(fit)
#?# Profiling...
#?# 2.5 % 97.5 %
#?# theta1 2.8448980 3.960982
#?# theta2 0.9433338 1.300629
#?# theta3 1.9147136 2.216113

#?# and 'fit' is:

#?# Call:
#?# mle(minuslogl = OU.lik, start = list(theta1 = 1, theta2 = 0.5,
#?# theta3 = 1), method = "L-BFGS-B", lower = c(-Inf, 0, 0))

#?# Coefficients:
#?# theta1 theta2 theta3
#?# 3.355322 1.106107 2.052815

#?# plus some graphic output

#?# Again, what's this telling me.

#?# This looks like a further example?
# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)




Please excuse the length of this email (and my lack of understanding)

Hope you can help and thanks.




Stephen Choularton Ph.D., FIoD


On 13/10/2010 2:41 AM, stefano iacus wrote:

just for completeness: OU process is gaussian and transitiion density is known in exact form. So maximum likelihood estimation works fine and I suggest to avoid GMM.

sde package contains exact transition density for this process (e.g. ?dcOU) which you can use to build the likelihood to pass to mle() function.

This example taken from the "inst" directory of the package sde. For the parametrization of the model see ?dcOU


# ex3.01.R
OU.lik <- function(theta1, theta2, theta3){
n <- length(X)
dt <- deltat(X)
-sum(dcOU(X[2:n], dt, X[1:(n-1)], c(theta1,theta2,theta3), log=TRUE))
}

require(stats4)
require(sde)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit
summary(fit)

# ex3.01.R (cont.)
prof <- profile(fit)
par(mfrow=c(1,3))
plot(prof)
par(mfrow=c(1,1))
vcov(fit)
confint(fit)

# ex3.01.R (cont.)
set.seed(123)
X <- sde.sim(model="OU", theta=c(3,1,2), N=1000, delta=1e-3)
mle(OU.lik, start=list(theta1=1, theta2=0.5, theta3=1),
method="L-BFGS-B", lower=c(-Inf,0,0)) -> fit2
summary(fit2)


I hope this helps out

stefano

On 12 Oct 2010, at 12:33, Bjorn Skogtro wrote:



Hi Stephen,

You could take a look at

http://sitmo.com/doc/Calibrating_the_Ornstein-Uhlenbeck_model

for the linear regression method, or take a look at the package "sde" which
contains some examples using GMM (not for the Ornstein-Uhlenbeck process,
though, only the CIR).

The half-life is given as log(2)/mean-reversion speed.

Do keep an eye on the partition of the time-axis, e.g. what frequency you
are using (daily, yearly) for interpreting the half-life.

BR,
Bjørn








------------------------------

Message: 2
Date: Tue, 12 Oct 2010 05:43:32 -0400
From: Sarbo
To: [hidden email]
Subject: Re: [R-SIG-Finance] Ornstein-Uhlenbeck
Message-ID:
Content-Type: text/plain; charset="utf-8"

By half-life, do you mean the speed of mean-reversion?

If so, there's a bit of algebraic tomfoolery that's required to
discretise the equation and then fit the data to it. I don't have the
time right now to go into all the details but it's not hard- you can
parameterise the process using simple linear regression. If you need
help with that I'll try and get back to you tonight about it.

On Tue, 2010-10-12 at 13:47 +1100, Stephen Choularton wrote:



Hi

Wonder if anyone could point me how I use this method to discover the
half life of a mean reverting process.

I am looking into pair trading and the time it takes for a
cointegrated pair to revert to the norm.

--
Stephen Choularton Ph.D., FIoD

9999 2226
0413 545 182


for insurance go to www.netinsure.com.au
for markets go to www.organicfoodmarkets.com.au


_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions


should go.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.html


-------------- next part --------------
A non-text attachment was scrubbed...
Name: CoS2010Winner.JPG
Type: image/jpeg
Size: 16091 bytes
Desc: not available
URL: <
https://stat.ethz.ch/pipermail/r-sig-finance/attachments/20101012/26e32fc7/attachment-0001.jpe


------------------------------

_______________________________________________
R-SIG-Finance mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-finance


End of R-SIG-Finance Digest, Vol 77, Issue 8
********************************************


[[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.



-----------------------------------
Stefano M. Iacus
Department of Economics,
Business and Statistics
University of Milan
Via Conservatorio, 7
I-20123 Milan - Italy
Ph.: +39 02 50321 461
Fax: +39 02 50321 505
http://www.economia.unimi.it/iacus
------------------------------------------------------------------------------------
Please don't send me Word or PowerPoint attachments if not
absolutely necessary. See:
http://www.gnu.org/philosophy/no-word-attachments.html

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com







_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.





No virus found in this incoming message.
Checked by AVG - www.avg.com




_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance --
Subscriber-posting only. If you want to post, subscribe first. -- Also
note that this is not the r-help list where general R questions should
go.
 		 	   		  
_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
No virus found in this incoming message. Checked by AVG - www.avg.com

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: Mean reversion

Yihao Lu aeolus_lu
In reply to this post by Sarbo

Thank you and thank Prof. Zivot

Best,
Yihao





From: [hidden email]
To: [hidden email]
Date: Tue, 19 Oct 2010 06:59:37 -0400
Subject: Re: [R-SIG-Finance] Mean reversion

Yihao,
 
Prof. Zivot is right. The ADF test isn't a great way to test for
mean-reversion; it's merely a way to test for stationarity to specify
the degrees in an ARIMA model.
 
I dug up some code from my misspent youth as a consultant which you
might find useful:
 
StochasticProcessTest <- function(PriceSeries, delta.t, p.value,
diagnostics = TRUE, from = NULL, to = NULL,
                                  by = NULL, currency = '($)'){
  require(stats)
  require(fSeries)
 
  m <- length(PriceSeries)
  S1 <- PriceSeries[1:(m-1)]
  S2 <- PriceSeries[2:m]
 
  #Fit a GBM process:
  Y <- diff(log(PriceSeries)); X <- S1
  GBMfit <- lm(Y ~ X)
 
  #Fit a mean-reverting GBM process:
  Y <- diff(PriceSeries) / S1; X <- log(S1)
  MRGBMfit <- lm(Y ~ X)
 
  #Fit a Vasicek (O-U) process:
  Y <- diff(PriceSeries); X <- S1
  Vfit <- lm(Y ~ X)
 
  #Fit a Cox-Ingersoll-Ross process:
  Y <- diff(PriceSeries) / sqrt(S1); X <- S1 / sqrt(S1)
  CIRfit <- lm(Y ~ X)
 
  #Gather the parameter estimates:
  kappa <- as.vector((c(0, -MRGBMfit$coef[2], -Vfit$coef[2], -CIRfit
$coef[2])) / delta.t)
  mu <- as.vector((c(mean(log(S2/S1)), MRGBMfit$coef[1] / kappa[2], Vfit
$coef[1] / kappa[3],
                     CIRfit$coef[1] / kappa[4])) / delta.t)
  sigma <- (c(sd(GBMfit$resid), sd(MRGBMfit$resid), sd(Vfit$resid),
sd(CIRfit$resid))) / sqrt(delta.t)
  tstat <- as.vector(c(summary(GBMfit)$coef[2,3],
summary(MRGBMfit)$coef[2,3], summary(Vfit)$coef[2,3],
                       summary(CIRfit)$coef[2,3]))
  fstat <- list(summary(GBMfit)$fstatistic,
summary(MRGBMfit)$fstatistic, summary(Vfit)$fstatistic,
                summary(CIRfit)$fstatistic)
  names(fstat) <- c('GBM', 'Mean.Revert.GBM', 'Vasicek', 'CIR')
  AICs <- c(AIC(GBMfit), AIC(MRGBMfit), AIC(Vfit), AIC(CIRfit))
  paramframe <- data.frame(rbind(kappa, mu, sigma, tstat, AICs))
  names(paramframe) <- names(fstat)
 
  #Now figure out what the actual process is, using AIC:
  crit <- ifelse(m > 30, qnorm(1 - p.value), qt(1 - p.value, df = n))
  tmp <- which.min(AICs)
  Processes <- names(fstat)
  Verdict <- Processes[tmp]
  FinalSummary <- switch(Verdict, GBM = summary(GBMfit), CIR =
summary(CIRfit), Mean.Revert.GBM = MRGBMfit,
                         Vasicek = Vfit)
  fitobj <- switch(Verdict, GBM = GBMfit, CIR = CIRfit, Mean.Revert.GBM
= MRGBMfit, Vasicek = Vfit)
  Output <- list(Parameters = paramframe, Critical.Value = crit, Verdict
= Verdict, FinalSummary = FinalSummary,
                 fstat = fstat, fitted.object = fitobj)
 
  if (diagnostics){
    op <- par(ask = TRUE)
    on.exit(op)
 
    if(!all(c(class(from), class(to)) == 'Date')){
      S <- timeSeries(PriceSeries)
    } else S = timeSeries(PriceSeries, seq(from, to, length.out = m))
 
    plot(S, type = 'l', xlab = 'Date', ylab = paste('Price Series',
currency), main = 'Time Series Plot of Data',
         lwd = 2, col = 'blue')
    rets <- returns(PriceSeries, 'continuous')[-1]
    hist(rets, xlab = paste('Log-Returns', currency), col = 'blue',
border = 'white', main = 'Histogram of Return Series',
         freq = FALSE)
    x <- seq(min(rets), max(rets), length.out = max(m, 1000))
    lines(x, dnorm(x, mean(rets), sd(rets)), col = 'magenta', lwd = 2)
    lines(density(rets), col = 'green', lwd = 2)
    legend('topright', legend = c('Log-Returns', 'Observed CDF',
'Gaussian Fit'), lwd = rep(2, 3),
                         col = c('blue', 'green', 'magenta'))
    plot(fitobj, 1:6)
  }
 
  return(Output)
 
}
   
 
(I don't claim that it's necessarily great code, but it does seem to
work.)
 
On Mon, 2010-10-18 at 21:35 -0400, Yihao Lu aeolus_lu wrote:
 
> I am doing rolling ADF test on some time series to check mean reversion. When I use short period rolling, I find the residue is not stationary at all. However, when I use horizon longer than 5 years, I find very significant stationary. On the other hand, I find the half life is only around 30 days.
> Is there anyone who can give me some possible explanation or guide me to some reference? thanks
>
> Best,
> Yihao

     
        [[alternative HTML version deleted]]

_______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-finance
-- Subscriber-posting only. If you want to post, subscribe first.
-- Also note that this is not the r-help list where general R questions should go.
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: cointegration

abhay
This post has NOT been accepted by the mailing list yet.
In reply to this post by Paul Teetor
Hi Paul,

I am very much impressed the way you explain your things. I have on question with regard to Mean Reversion I got various pairs my idea is to identify the indictor which will give me an indication of reversion. Can you highlight any method which would probably help me to identify such reversion.

Thanks
Loading...