R/ar1est.R

ar1est <- function(z, method=c("MLE", "LSE") ){
    n <- length(z)
    stopifnot(n > 3)
    WhichMethod <- match.arg(method)
    m <- mean(z)
    MLEQ <- identical(WhichMethod, "MLE")
    if (MLEQ) {
        mz <- rep(m, n)
        a <- sum((z - mz)^2)
        b <- sum((z[-1] - mz[-1]) * (z[-n] - mz[-n]))
        c <- sum((z[c(-1, -n)] - mz[c(-1, -n)])^2)
        i <- complex(1, 0, 1)
        x <- ((-16) * b^3 + 18 * a * b * c + 24 * b^3 * n - 27 * 
            a * b * c * n - 9 * b * c^2 * n - 12 * b^3 * n^2 + 9 * 
            a * b * c * n^2 + 27 * b * c^2 * n^2 + 2 * b^3 * n^3 - 
            18 * b * c^2 * n^3)
        y <- (-b^2 * (-2 + n)^2 + 3 * c * (-1 + n) * (-a - c * n))
        f <- complex(1, x^2 + 4 * y^3, 0)
        z <- (x + sqrt(f))^(1/3)
        g <- x^2 + 4 * y^3
        z1 <- (x + (-g)^(1/2) * i)^(1/3)
        part1 <- (n - 2) * b/(3 * c * (n - 1))
        part2 <- (1 - sqrt(3) * i) * y/(3 * 2^(2/3) * c * (n - 1) * z)
        part3 <- (1 + sqrt(3) * i) * z/(6 * 2^(1/3) * c * (n - 1))
        rho <- Re(part1 + part2 - part3)
     }
    else  {
        m1 <- mean(z[-1])
        m2 <- mean(z[-n])
        rho <- sum((z[-1]-m1)*(z[-n]-m2))/sum((z[-n]-m2)*(z[-n]-m2))
    }
    rho
}

Try the mleur package in your browser

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

mleur documentation built on May 1, 2019, 7:31 p.m.