inst/scripts/tseries/AR1Gumbel.R

In ROptRegTS: Optimally Robust Estimation for Regression-Type Models

###############################################################################
# AR(1) with Gumbel innovations -- *=c, t=alpha=1
###############################################################################

require(distrEx)

###############################################################################
# Computation of centering constant
###############################################################################
.AR1getah <- function(a, c0){
    c0 <- c0^2/(1-c0^2)
    c1 <- max(a-c0,0)
    c2 <- c0 + a

    return((1 + c1 - (a-c0))*exp(-c1) - exp(-c2) - c0)
}

.AR1geta <- function(c0, eps){
    uniroot(.AR1getah, lower = 0, upper = 1.01, tol = eps, c0 = c0)$root 
}


###############################################################################
# Computation of optimal clipping bound
###############################################################################
.AR1getc0 <- function(c0, H, rad, FFT, eps){
    integrandc0 <- function(x, c0){
        c.H <- c0/pmax(1e-10,abs(x))
        c.H1 <- sqrt(c.H/(1+c.H))        
        if(FFT)
            a.H <- sapply(c.H1, .AR1geta, eps = eps)
        else
            a.H <- .AR1geta.spline(c.H1)

        c1 <- pmax(a.H-c.H,0)
        c2 <- a.H + c.H
        r.H <- abs(x)*pmax(((1 + c1 - (a.H-c.H))*exp(-c1) + exp(-c2) + (a.H-c.H) - 1), 0)
    }
    
    if(FFT)
        r1 <- E(H, integrandc0, c0 = c0)
    else
        r1 <- mean(integrandc0(x = H, c0 = c0))
    r <- sqrt(r1/c0)
    
    return(rad-r)
}

###############################################################################
# Computation of standardizing constant
###############################################################################
.AR1getA <- function(x, c0, FFT, eps){
    c.H <- c0/pmax(1e-10,abs(x))
    c.H1 <- sqrt(c.H/(1+c.H))
    if(FFT)
       a.H <- sapply(c.H1, .AR1geta, eps = eps)
    else
       a.H <- .AR1geta.spline(c.H1)
    c1 <- pmax(a.H-c.H,0)
    c2 <- a.H + c.H
    aa.H <- (c1*(2 + c1 - (a.H-c.H)) + 2 - (a.H-c.H))*exp(-c1) - (2+c2)*exp(-c2) - c.H

    return(x^2*aa.H)
}


###############################################################################
# computation of optimal IC
###############################################################################
AR1Gumbel <- function(rad, phi, shift = 0, FFT = TRUE, delta = .Machine$double.eps^0.5, 
                     simn = 100000, offs = 10000, x.start = 0, 
                     eps = .Machine$double.eps^0.5, upper = 1000){  
    if(length(rad) != 1)
        stop("'rad' has to be of length 1")
    if(rad <= 0)
        stop("'rad' has to be in (0,Inf)")
    if(length(phi) != 1)
        stop("'phi' has to be of length 1")
    if((phi <= -1)||(phi >= 1))
        stop("'phi' has to be in (-1,1)")
    if(FFT){
        Innov <- Gumbel(loc = digamma(1))
        convpow <- ceiling(log(delta)/log(abs(phi)))
        X <- Innov
        for(i in 1:convpow){
            X.alt <- X
            X <- shift + X.alt + (-phi)^i*Innov
        }
        rm(X.alt)
        H <- -X
        rm(X)
    }else{
        x <- x.start
        X <- numeric(simn)
        for(i in 1:(simn + offs)){
            x <- shift - phi*x + rgumbel(1, loc = digamma(1))
            if(i > offs) X[i-offs] <- x
        }
        H <- -X
        rm(X)
        if(!exists(".AR1geta.spline")){
            c.vct <- seq(from = 1e-6, to = 0.99, length = 10000)
            a.vct <- sapply(c.vct, .AR1geta, eps = eps)
            assign(".AR1geta.spline", splinefun(c.vct, a.vct), env = .GlobalEnv)
        }
    }

    c0 <- uniroot(.AR1getc0, lower = 1e-4, upper = upper, tol = eps, H = H, 
                  rad = rad, eps = eps, FFT = FFT)$root

    if(FFT)
        A <- 1/E(H, .AR1getA, c0 = c0, FFT = FFT, eps = eps)
    else{
        A <- 1/mean(.AR1getA(x = H, c0 = c0, FFT = FFT, eps = eps))
        old <- distroptions("DistrResolution")
        distroptions("DistrResolution", 1e-12)
        H <- DiscreteDistribution(H)
        distroptions("DistrResolution", old)
    }
    b <- A*c0

    return(list(A = A, b = b, H = H))
}