Regression technique with correlated errors, based on Cochrane-Orcutt (http://en.wikipedia.org/wiki/Cochrane%E2%80%93Orcutt_estimation)

tolleyshuffle.R

## Regression with correlated errors

tolley.shuffle <- function(model,k=1e-5){
    if(class(model)!="lm")
        stop("Object 'model' must be of class 'lm'.")

    uncorrelate <- function(model) {
        r <- residuals(model)
        n <- length(r)
        y <- model$model[,1]
        z <- model$model[,2:ncol(model$model)]
    
        ## Determine p and q in the ARMA(p,q)
        ACF <- acf(r)$acf[1:10]
        PACF <- pacf(r)$acf[1:10]
    
        critical <- 2/sqrt(n)
        a <- sum(ACF > critical)
        p <- sum(PACF > critical)
    
        if (a >= 4 & p <= 2)
            order <- c(p,0,0)
        else if (a <= 2 & p >= 4)
            order <- c(0,0,a)
        else order <- c(1,0,1)
    
        ## Fit the ARMA model and extract AR and MA components
        mod.r <- arima(r,order=order)    
        phi <- mod.r$coef[substr(names(mod.r$coef),1,1)=='a']
        theta <- mod.r$coef[substr(names(mod.r$coef),1,1)=='m']
    
        ## Transform functions
        transform <- function(x,param) {
            if (sum(param) > 0) {
                x.t <- matrix(0,length(param),length(x))
                for(i in 1:length(param)) {
                    x.t[i,] <- lag(ts(x),-i)*param[i]
                }
                x.t <- colSums(x.t)
            } else x.t <- x
            return(x.t)
        }
        
        y.phi <- transform(y,phi)
        y.theta <- transform(y,theta)
        z.phi <- transform(z,phi)
        z.theta <- transform(z,theta)
    
        if (all(y.phi==y))
            new.y <- y.theta
        else if (all(y.theta==y))
            new.y <- y.phi
        else new.y <- y.phi/y.theta
    
        if (all(z.phi==z))
            new.z <- z.theta
        else if (all(z.theta==z))
            new.z <- z.phi
        else new.z <- z.phi/z.theta
    
        new.mod <- lm(new.y~new.z)
        return(new.mod)
    }

    ## Sets everything up for iterations until convergence
    mod <- model
    par.sums <- NULL
    difference <- j <- 1
    while(difference > k) {
        j <- j+1
        mod <- uncorrelate(mod)
        par.sums <- c(par.sums,sum(mod$coef))
        difference <- par.sums[j]-par.sums[j-1]
    }
    return(list(mod,j))
}