R/fdSperio.R

Defines functions fdSperio

Documented in fdSperio

#### by Valderio Reisen  -- Dec.2005--
#### Tweaks by MM

## MM(FIXME): This is "in parallel"  to fdGPH() , see ./fdGPH.R

fdSperio <- function(x, bandw.exp = 0.5, beta = 0.9)
{
    if(NCOL(x) > 1) stop("only implemented for univariate time series")
    x <- as.numeric(na.fail(as.ts(x)))
    if (any(is.na(x))) stop("NAs in x")
    n <- length(x)
    ## Compute "smoothed" periodogram -- MM (FIXME): use  spec.pgram() !
    g <- trunc(n^bandw.exp)
    j <- 1:g
    kk <- 1:(n-1)
    w <- 2*pi*j/n
    x <- x - mean(x)
    var.x <- sum(x^2)/n # not /(n-1)
    cov.x <- numeric(n-1)
    for (k in kk)
        cov.x[k] <- sum(x[1:(n-k)] * x[(1+k):n]) / n

    M <- trunc(n^beta)
    M2 <- M %/% 2
    pw <- numeric(n-1)
    for (k in kk) {
        A_k <- k/M
        pw[k] <-
            if (k <= M2)  1 - 6*A_k^2 *(1 - A_k)
            else if (k <= M) 2*(1 - A_k)^3 else 0
    }
    periodogram <- numeric(g)
    for (i in 1:g) # unscaled (will scale below)
        periodogram[i] <- var.x + 2*sum(cov.x* pw * cos(w[i]*kk))

    pos <- j[periodogram > 0]
    y.reg <- log(periodogram[pos] / (2*pi))
    x.reg <- 2*log(2*sin(w[pos]/2)) ## = log( (2*sin(..)) ^ 2)
    fit <- lm.fit(cbind(1, x.reg), y.reg)
    d.GPH <- coef(fit)[["x.reg"]]
    x.r2 <- sum((x.reg - mean(x.reg))^2)
    var.d <- (0.539285*M/n)/ x.r2
    var.reg <- sum(resid(fit)^2) / ((g - 1) * x.r2)
    list(d = -d.GPH, sd.as = sqrt(var.d), sd.reg = sqrt(var.reg))
}

Try the fracdiff package in your browser

Any scripts or data that you put into this service are public.

fracdiff documentation built on Nov. 1, 2022, 1:06 a.m.