haarhist <- function(dat, j = 2)
{
  xgrid <- (0:10000)/2^j
  xgrid <- xgrid[xgrid >= min(dat)]
  xgrid <- xgrid[xgrid <= max(dat)]
  xgrid <- c(xgrid[1] - 2^( - j), xgrid, xgrid[length(xgrid)] + 2^( - j))
  cjk <- rep(0, length(xgrid) - 1)
  for(i in 1:length(cjk))
    cjk[i] <- ((sum(dat <= xgrid[i + 1]) - sum(dat <= xgrid[i])) * 
	       2^j)/length(dat)
  plot(rep(xgrid, rep(2, length(xgrid))), c(0, rep(cjk, rep(2, length(cjk))), 
	       0), type = "l")
  mtext(side = 3, line = 0.1, paste("j = ", j))
  NULL
}