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 <[EMAIL PROTECTED]>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 <[EMAIL PROTECTED]>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 <[EMAIL PROTECTED]>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@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.
