Theil-Sen regression

tsreg.R

# Ported by David Holstius <dholstius@baaqmd.gov>
# from http://skip.ucsc.edu/leslie_MOUSE/programs/plotting/tsreg.r

# Module tsreg
# Author:  E. A. Houseman
# Last update July 2004

# AR(q) time series regression assuming regular intervals
# Support for cholesky residuals [Houseman, Ryan, Coull (2004, JASA)]

tsreg <- function (model, data, ord = 1, maxiter = 10, tol = sqrt(.Machine$double.eps)){
  
  y <- model.extract(model.frame(model,data),"response") # Get response
  X <- model.matrix(model,data) # Get covariate matrix
  
  dd <- dim(X)
  n <- dd[1]
  p <- dd[2]
  
  beta <- solve(t(X) %*% X, t(X) %*% y)
  e <- y - as.vector(X %*% beta)
  
  if(ord == 0){
    
    sigma <- sqrt( sum(e*e)/(n-p) )
    V.beta <- sigma * solve(t(X) %*% X)
    V.Sigma <- matrix(2 * sigma^4 * n/(n-p)^2, 1, 1)
    o <- list(beta = beta, rho = numeric(0), iters = 0, 
              ord = ord, V.beta = V.beta, V.Sigma = V.Sigma, 
              R.elems  =  NULL, dR.elems  =  NULL,
              sigma = sigma, y = y, X = X, shifts = NULL, call = match.call())
    class(o) <- "tsreg"
    return(o)
  }
  
  shifts <- tsreg.makeShifts(n)  
  Z <- tsreg.makeZ(e,ord,n)
  rho0 <- rho <- solve(t(Z) %*% Z, t(Z) %*% e[-(1:ord)])
  
  for(i in 1:maxiter){
    R.elems <- tsreg.makeCorrElems(rho, shifts)
    R <- tsreg.makeCorr(R.elems)
    R.qr <- qr(R)
    Xt <- t(solve(R.qr,X))
    XtX <- Xt %*% X
    beta <- solve(XtX, Xt %*% y)
    e <- y - as.vector(X %*% beta)
    Z <- tsreg.makeZ(e,ord,n)
    
    rho0 <- rho
    ZtZ <- t(Z) %*% Z
    rho <- solve(ZtZ, t(Z) %*% e[-(1:ord)])  
    
    if(max(abs(rho-rho0))<tol) break;
  }
  
  sigma <- as.vector(sqrt(e %*% solve(R.qr,e)/n))
  sigma2 <- sigma*sigma
  
  bnames <- rownames(beta)
  beta <- as.vector(beta)
  names(beta) <- bnames
  
  V.beta <- sigma2*solve(XtX)
  
  dR.elems <- tsreg.makeDCorrElems(rho, shifts)
  RiR. <- list()
  Info <- matrix(0, ord+1, ord+1)
  
  for(j in 1:ord){
    R. <- tsreg.makeCorr(dR.elems[[j]])
    diag(R.) <- 0
    RiR.[[j]] <- solve(R.qr, R.)
    
    for(k in 1:j){
      if(j == k) Info[j,j] <- sum(diag( RiR.[[j]] %*% RiR.[[j]] ))
      else Info[j,k] <- Info[k,j] <- sum(diag(RiR.[[k]] %*% RiR.[[j]]))
    }
    Info[j,ord+1] <- Info[ord+1,j] <- sum(diag(RiR.[[j]])) / sigma2
  }
  Info[ord+1,ord+1] <- n / sigma2 / sigma2
  
  V.Sigma <- 2 * solve(Info)
  
  obj <- list(beta = beta, rho = as.vector(rho), iters = i, 
            ord = ord, V.beta = V.beta, V.Sigma = V.Sigma, 
            R.elems  =  R.elems, dR.elems  =  dR.elems, 
            sigma = sigma, y = y, X = X, shifts = shifts, call = match.call())
  class(obj) <- "tsreg"
  return(obj)
}

# Redefine chol as a method
chol <- function (x, ...) UseMethod("chol")
chol.matrix <- base::chol

# Summary for tsreg
summary.tsreg <- function (o, alpha = 0.05){
  beta <- cbind(est = o$beta, SE = sqrt(diag(o$V.beta)))
  Sigma <- cbind(est = c(o$rho,o$sigma^2), SE = sqrt(diag(o$V.Sigma)))
  
  crit <- abs(qnorm(alpha/2))
  beta <- cbind(beta, beta[,1] - crit*beta[,2])
  beta <- cbind(beta, beta[,1] + crit*beta[,2])
  Sigma <- cbind(Sigma, Sigma[,1] - crit*Sigma[,2])
  Sigma <- cbind(Sigma, Sigma[,1] + crit*Sigma[,2])
  sig.parm <- dim(Sigma)[1]
  
  # Base CI's for sigma^2 on log transformation
  Sigma[sig.parm,3] <- exp(log(Sigma[sig.parm,1])-crit*Sigma[sig.parm,2]/Sigma[sig.parm,1])
  Sigma[sig.parm,4] <- exp(log(Sigma[sig.parm,1])+crit*Sigma[sig.parm,2]/Sigma[sig.parm,1])
  
  colnames(beta)[3:4] <- paste((1-alpha)*100,"% CI ",c("lo","hi"),sep = "")
  colnames(Sigma)[3:4] <- paste((1-alpha)*100,"% CI ",c("lo","hi"),sep = "")
  
  if(o$ord == 0) rownames(Sigma) <- "sigma^2"
  else rownames(Sigma) <- c(paste("rho[",1:length(o$rho),"]",sep = ""), "sigma^2")
  
  so <- list(call = o$call, beta = beta, Sigma = Sigma)
  class(so) <- "tsregsummary"
  so
}

print.tsreg <- function (fit){
  cat("Output from 'tsreg'\n")
  print(fit$call)
  cat("\nIterations:",fit$iter,"\n") 
  cat("\nRegression parameters:\n")
  print(fit$beta)
  cat("\nCorrelation parameters:\n")
  print(fit$rho)
  cat("\nCovariance scale (sigma^2)\n") 
  print(fit$sigma^2)
}

print.tsregsummary <- function (sfit){
  cat("Summary of output from 'tsreg'\n")
  print(sfit$call)
  cat("\nRegression parameters:\n")
  print(sfit$beta)
  cat("\nCovariance parameters:\n")
  print(sfit$Sigma)
  cat("\n")
}

# Cholesky residuals for tsreg solution
chol.tsreg <- function (o){
  if(o$ord == 0) R <- diag(length(o$y))
  else R <- tsreg.makeCorr(o$R.elems)
  
  R.qr <- qr(R)
  R.inv <- solve(R.qr)
  
  Lt <- chol(R.inv) / o$sigma
  z <- as.vector(Lt %*% (o$y - o$X %*% o$beta))
  
  ord <- o$ord
  p <- length(o$beta)
  n <- length(z)
  
  ds2 <- matrix(0,n,ord+1)
  
  ds2[,ord+1] <- -diag(R.inv)/o$sigma^2
  if(ord>0) for(r in 1:ord){
    R. <- tsreg.makeCorr(o$dR.elems[[r]])
    diag(R.) <- 0
    ds2[,r] <- - diag(Lt %*% R. %*% t(Lt))
  }
  ds2 <- ds2/o$sigma^2
  dm <- Lt %*% o$X
  aux <- new.env()
  assign("dm", dm, env = aux)
  assign("ds2", ds2, env = aux)
  assign("V.beta", o$V.beta, env = aux)
  assign("V.Sigma", o$V.Sigma, env = aux)
  attr(z,"aux") <- aux
  class(z) <- "tsregchol" 
  z
}

# Summary of Cholesky residuals for getting confidence bands
summary.tsregchol <- function (z, x = z, alpha = 0.05){
  p <- pnorm(x)
  n <- length(z)
  
  dm <- apply(get("dm", env = attr(z,"aux")), 2, mean)
  ds2 <- apply(get("ds2", env = attr(z,"aux")), 2, mean)
  V.beta <- get("V.beta", env = attr(z,"aux"))
  V.Sigma <- get("V.Sigma", env = attr(z,"aux"))
  
  delta.m <- t(dm) %*% V.beta %*% dm
  delta.s2 <- t(ds2) %*% V.Sigma %*% ds2
  
  dnormx <- dnorm(x)
  var.correction <- as.vector(dnormx*dnormx*delta.m + 
                                x*x*dnormx*dnormx*delta.s2/4)
  
  var.p <- p*(1-p)/n 
  se <- sqrt(var.p - var.correction)
  z.alpha <- qnorm(alpha/2)
  ci.lo <- p - z.alpha * se
  ci.hi <- p + z.alpha * se
  
  sort <- order(x)
  cbind(p,se,ci.lo,ci.hi,zci.lo = qnorm(ci.lo),zci.hi = qnorm(ci.hi),index = sort)[sort,]
}

# QQ plot with confidence "bands"
qqnorm.tsregchol <- function (z,xlab = NULL,ylab = NULL,col = 2,...){
  band <- suppressWarnings( summary(z) )
  qqnorm(as.vector(z),xlab = xlab,ylab = ylab)
  n <- length(z)
  qz <- qnorm(band[,1])
  lines(band[,5],qz,col = col,...)
  lines(band[,6],qz,col = col,...)
}

qqnorm.tsreg <- function (fit,xlab = NULL,ylab = NULL,col = 2,...){
  qqnorm(chol(fit),xlab = xlab,ylab = ylab,col = col,...)
}


#########################
#### SUPPORTING FUNCTIONS
#########################

# Make the covariate matrix for LS regression to get rho parameters
# [see Durbin (1960)]
tsreg.makeZ <- function (e, ord, n){
  Z <- matrix(0,n-ord,ord)
  for(j in 1:ord){
    Z[,j] <- e[(ord-j+1):(n-j)]
  }
  Z   
}

# Get shift matrices for constructing correlations
tsreg.makeShifts <- function (n){
  E <- list()
  E[[1]] <- diag(n)
  for(i in 1:n) {
    e <- rep(0,n)
    if(i>1) e[i-1] <- 1
    E[[i+1]] <- rbind(e,E[[i]][-n,])
  }
  E
}

# Make elements of a correlation matrix
tsreg.makeCorrElems <- function (rho, shifts){
  shiftmat <- tsreg.makeShiftMat(rho, shifts) 
  solve(shiftmat$E0, shiftmat$Rho)
}

# Sum of shift matrices for constructing correlation matrix
tsreg.makeShiftMat <- function (rho, shifts){
  E0 <- shifts[[1]]
  n <- length(shifts)-1
  Rho <- c(rho, rep(0,n-length(rho)))
  for(i in 1:n){
    E0 <- E0 - Rho[i]*shifts[[i+1]]
  }
  list(E0 = E0, Rho = Rho)
}

# Make a correlation matrix
tsreg.makeCorr <- function (elems){
  n <- length(elems)
  R <- matrix(0,n,n)
  for(i in 1:(n-1)){
    R[i,i+(1:(n-i))] <- elems[1:(n-i)]
  }
  R <- R + t(R)
  diag(R) <- 1
  R
}

# Make elements of correlation derivatives
tsreg.makeDCorrElems <- function (rho, shifts){
  ord <- length(rho)
  dR.elems <- list()
  shiftmat <- tsreg.makeShiftMat(rho, shifts)
  E0.qr <- qr(shiftmat$E0)
  for(r in 1:ord){
    E.dot <- shifts[[r+1]]
    dR.elems[[r]] <- as.vector(solve(E0.qr)[,r] + 
                                 solve(E0.qr, E.dot) %*% solve(E0.qr, shiftmat$Rho))
  }
  dR.elems
}