R/DLResiduals.R

In ltsa: Linear Time Series Analysis

`DLResiduals` <-
function(r,z,useC=TRUE, StandardizedQ=TRUE)
{
    n <- length(z)
    if (n!=length(r)) stop("arguments have unequal length")
    EPS <- .Machine$double.eps # 1+EPS==1, machine epsilon
    if (StandardizedQ)
        stanQ <- 1
    else
        stanQ <- 0
    if (useC){
      out<-.C("DLResid", as.double(z), res=as.double(numeric(length(z))), as.integer(length(z)), as.double(r), 
            LogL=as.double(1), EPS, as.integer(stanQ), fault = as.integer(1),
      PACKAGE="ltsa" )
    fault<-out$fault
    if (fault == 1) 
            stop("error: sequence not p.d.")
    res<-out$res
    }
    else {
        error <- numeric(n)
        sigmasq <- numeric(n)
        error[1] <- z[1]
        sigmasq[1] <- r[1]
        phi <- r[2]/r[1]
        error[2] <- z[2] - phi * z[1]
        sigmasqkm1 <- r[1] * (1 - phi^2)
        sigmasq[2] <- sigmasqkm1
        for(k in 2:(n - 1)) {
            if (sigmasqkm1<EPS) stop("r is not a p.d. sequence")
            phikk <- (r[k + 1] - phi %*% rev(r[2:k]))/sigmasqkm1
            sigmasqk <- sigmasqkm1 * (1 - phikk^2)
            phinew <- phi - phikk * rev(phi)
            phi <- c(phinew, phikk)
            sigmasqkm1 <- sigmasqk
            error[k + 1] <- z[k + 1] - crossprod(phi, rev(z[1:k]))
            sigmasq[k + 1] <- sigmasqk
            }
        if (StandardizedQ)
            res<-error/sqrt(sigmasq)
        else
            res<-error
        }
    res
}