Re: [R] S-PLUS code in R

From: Fotis Papailias <fpapailias_at_gmail.com>
Date: Sat, 26 Jul 2008 13:52:55 +0100

Dear R-Users,

the problem I was facing it's solved. Actually, the code was ok from the first time.

You define it as:

lw <- function(h, x, im)
{

    peri1 <- per(x)
    len <- length(x)
    m <- len/im
    peri <- peri1[2:(m+1)]
    z <- c(1:m)
    freq <- ((2*pi)/len) * z
    result <- log(sum(freq^(2*h-1)*peri))-(2*d)/m * sum(log(freq)) }

and after you have input (or simulate) your series (say x) you optimize the function as

k <- optimize(lw, interval=c(0.5, 1), x=x1, im=2)

and you get the estimated value you are interested in as *
*k2 <- k$minimum

Just for the record, so it this message can be useful for future mailing list searches, the above function is the so called Local Whittle or Gaussian Semiparametric Estimation first introduced by Robinson (1995) and studied by Taqqu (1995-...) who is the author of the original S-PLUS code of my above function.

You can use it when you want to estimate the fractional exponent (or long memory parameter) in time series exhibit long memory. However you have to notice that, this function takes in account the Hurst exponent (denoted as "H" or "h" in the code) and the long memory exponent "d", consequently, when you want to estimated an ARFIMA (p, *d*, q) series where d=0.2, the "e2" value that the code will return you is the estimated value of *H* which is going to be around 0.7 . The relation between H and d is defined as H=d+0.5, hence when you get "e2" the one that you are really looking for is

e.lwhittle <- e2-0.5

Take care when you are changing the "im" parameter and visit the page of Taqqu as mentioned somewhere below to find more theory and proofs about the local whittle and about what all the parameters in function lw are doing.

Best wishes,

Fotis

On Sat, Jul 26, 2008 at 1:28 PM, Fotis Papailias <fpapailias_at_gmail.com>wrote:

> I made some experiments and I conclude that if I change the interval (e.g.
> say c(-100,100)!!!) then the minimum value changes in respect to the
> simulated data.
>
> So, probably the problem lies in the "optimal" values of the interval?
>
>
>
> On Sat, Jul 26, 2008 at 1:00 PM, Fotis Papailias <fpapailias_at_gmail.com>wrote:
>
>> Hi,
>>
>> thanks for your message.
>>
>> You mean to rewrite the function like that:
>>
>> lw <- function(d, x, im)
>> {
>> peri1 <- per(x)
>> len <- length(x)
>> m <- len/im
>> peri <- peri1[2:(m+1)]
>> z <- c(1:m)
>> freq <- ((2*pi)/len) * z
>> result <- log(sum(freq^(2*d-1)*peri))-(2*d)/m * sum(log(freq))
>> }
>>
>>
>> and then to use:
>>
>> target <- function(d) lw(x, d, im=2)
>> k <- optimize(target, interval=c(0, 0.5))
>>
>> or just,
>>
>> k <- optimize(lw, x, im)
>>
>> But I still get the same values
>>
>> $minimum
>> [1] 0.4999542
>>
>> $objective
>> [1] 2.739509
>>
>> after the optimization, no matter my "x" series. It seems that I am still
>> doing something wrong in the optimize() or in the result() function inside
>> the loop. But comparing it to the original S-PLUS code, it seems like it is
>> fine.
>>
>> Send me if you have any other ideas please.
>>
>> Thanks again.
>>
>> fotis
>>
>>
>> On Sat, Jul 26, 2008 at 12:48 PM, Duncan Murdoch <murdoch_at_stats.uwo.ca>wrote:
>>
>>> On 26/07/2008 7:40 AM, Fotis Papailias wrote:
>>>
>>>> Dear R-users,
>>>>
>>>> I have sent another mail some hour ago about a matlab Code I was trying
>>>> to
>>>> translate in R.
>>>>
>>>> Actually I have found a simpler code originally written in S-PLUS for
>>>> the
>>>> same function.
>>>> Author's page -> http://math.bu.edu/people/murad/methods/locwhitt/
>>>>
>>>> =============================================================
>>>>
>>>> rfunc_function(h, len, im, peri)
>>>> # h -- Starting H value for minimization.
>>>> # len -- Length of time series.
>>>> # im -- Use only len/im frequencies.
>>>> # peri -- Periodogram of data.
>>>> {
>>>> m <- len %/% im
>>>> peri <- peri[2:(m + 1)]
>>>> z <- c(1:m)
>>>> freq <- (2 * pi)/len * z
>>>> result <- log(sum(freq^(2 * h - 1) * peri)) - (2 * h)/m *
>>>> sum(log(freq)
>>>> ) # cat("H = ", h, "R = ", result, "\n")
>>>> drop(result)
>>>> }
>>>>
>>>>
>>>> locwhitt_function(data, h = 0.5, im = 2)
>>>> # data -- Time series.
>>>> # h -- Starting H value for minimization.
>>>> # im -- Use only N/im frequencies where N is length of series.
>>>>
>>>> {
>>>> peri <- per(data)
>>>> len <- length(data)
>>>> return(nlminb(start = h, obj = rfunc, len = len, im = im, peri =
>>>> peri)$
>>>> parameters)
>>>> }
>>>> ===============================================================
>>>>
>>>> The author who has written the above S-PLUS code uses two functions
>>>> (with
>>>> the locwhitt_function he lets the user to define the data and the
>>>> parameters
>>>> and with the rfunc_function he does the minimization.)
>>>>
>>>> Mine translation is in R is:
>>>>
>>>> where I use a joint function compared to the the above author
>>>>
>>>>
>>>> ================================================================
>>>>
>>>> lw <- function(x, d, im)
>>>> {
>>>> peri1 <- per(x)
>>>> len <- length(x)
>>>> m <- len/im
>>>> peri <- peri1[2:(m+1)]
>>>> z <- c(1:m)
>>>> freq <- ((2*pi)/len) * z
>>>> result <- log(sum(freq^(2*d-1)*peri))-(2*d)/m * sum(log(freq))
>>>> }
>>>>
>>>> =================================================================
>>>>
>>>> which seems to run ok.
>>>>
>>>> But when I do
>>>>
>>>> k <- optimize(lw, x, im=2, interval=c(0, 0.5))
>>>>
>>>> I always get the same result no matter the (simulated) data in x!
>>>>
>>>> The parameter of interest to be minimized is "d". Does anyone know how
>>>> to
>>>> edit the function "optimize" so it can work properly?
>>>>
>>>
>>> optimize() is fine, but the way you're calling it is not. It optimizes a
>>> function over the first argument. So you could rewrite lw to put d first,
>>> or write a new function which calls it, e.g.
>>>
>>> target <- function(d) lw(x, d, im)
>>>
>>> and then
>>>
>>> optimize(target, interval=c(0, 0.5))
>>>
>>> Because target is defined in the global environment, it will look there
>>> for x and im, and you don't need to pass them as arguments: unless x and im
>>> aren't defined there too!
>>>
>>> Duncan Murdoch
>>>
>>
>>
>>
>> --
>> fp
>>
>
>
>
> --
> fp
>

-- 
fp

	[[alternative HTML version deleted]]

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Received on Sat 26 Jul 2008 - 12:56:43 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sat 26 Jul 2008 - 13:32:22 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive