###################

## Preliminaries ##

###################

## packages

library("coin")

library("strucchange")

library("lattice")

## random seed

rseed <- 20061103

## theme for lattice/trellis graphics

trellis.par.set(theme = canonical.theme(color = FALSE))

######################

## Boston homicides ##

######################

## time series plot

data("BostonHomicide", package = "strucchange")

hom_month <- zoo(as.vector(BostonHomicide$homicides), as.yearmon(as.vector(time(BostonHomicide$homicides))))

hom_year <- aggregate(hom_month, function(x) as.numeric(floor(x)), mean)

plot(hom_month, col = grey(0.7), lwd = 2, ylab = "Number of homicides", xlab = "Time")

lines(hom_year, type = "b")

## monthly data

hom_month <- data.frame(

  log_homicides = log(coredata(hom_month) + 0.5),

  time = as.numeric(time(hom_month)))

## asymptotic unconditional test

sctest(gefp(log_homicides ~ 1, data = hom_month), functional = supLM(0.1))

## asymptotic conditional test

set.seed(rseed)

maxstat_test(log_homicides ~ time, data = hom_month, minprob = 0.1)

## approximate conditional test

set.seed(rseed)

maxstat_test(log_homicides ~ time, data = hom_month, minprob = 0.1, distribution = approximate(9999))

## annual data

hom_year <- data.frame(

  homicides = coredata(hom_year),

  time = time(hom_year))

## asymptotic unconditional test

sctest(gefp(homicides ~ 1, data = hom_year), functional = supLM(0.1))

## asymptotic conditional test

set.seed(rseed)

maxstat_test(homicides ~ time, data = hom_year, minprob = 0.1)

## approximate conditional test

set.seed(rseed)

maxstat_test(homicides ~ time, data = hom_year, minprob = 0.1, distribution = approximate(9999))

## note that the manuscript computes the exact conditional p-value

## by exhaustive search (rather than approximating via 10,000 simulations)

###########################

## Hiring discrimination ##

###########################

## data

hire <- cbind(c(2, 0, 0, 0, 5, 14), c(427, 86, 104, 180, 111, 59))

hire <- data.frame(

  resp = factor(unlist(sapply(1:nrow(hire), function(i) rep(c("female", "male"), hire[i,])))),

  time = as.numeric(rep(1991:1996, rowSums(hire))))

## visualization

set.seed(rseed)

maxstat_test(resp ~ time, data = hire, distribution = approximate(9999), minprob = 0.05)

################

## CO2 reflux ##

################

## data

data("CWD", package = "coin")

lwood <- reshape(CWD, varying = paste("sample", c(2:4,6:8), sep = ""),

  timevar = "tree", direction = "long", v.names =  "sample")

lwood$tree <- factor(lwood$tree)

## visualization

print(xyplot(sample ~ time | tree, data = lwood, type = "b",

  scales = list(y = list(relation = "free")),

  layout = c(3, 2), xlab = "Time", ylab = expression(paste(CO[2], " reflux")),

  as.table = TRUE))

## test

(cwd_mt <- maxstat_test(sample2 + sample3 + sample4 + sample6 + sample7 + sample8 ~ trend, 

  data = CWD, distribution = approximate(1e5)))

## maximally selected statistics with 5% critical value

cwd_st <- statistic(cwd_mt, "standardized")

cwd_st <- data.frame(cwd_st, time = CWD$time[1] + CWD$trend[2:11])

cwd_st <- reshape(cwd_st, varying = paste("sample", c(2:4,6:8), sep = ""),

  timevar = "tree", direction = "long", v.names = "sample")

cwd_st$tree <- factor(cwd_st$tree)

cwd_q <- qperm(cwd_mt, 0.95)

print(xyplot(sample ~ time | tree, data = cwd_st, type = "b",

  panel = function(...) {

    panel.xyplot(...)

    panel.abline(h = c(-cwd_q, cwd_q), col = 2)

    panel.abline(h = 0, col = "gray")

  },

  layout = c(3, 2), xlab = "Time", ylab = "Test statistics", ylim = c(-3.5, 3.5),

  as.table = TRUE))

##################################

## Dow Jones Industrial Average ##

##################################

## data

data("DJIA", package = "strucchange")

djia <- diff(log(DJIA)) * 100

djia_res <- coredata(djia - mean(djia))

djia_trafo <- cbind(Intercept = djia_res, Variance = djia_res^2 - mean(djia_res^2))

djia_time <- time(djia)

## visualization

plot(djia, xlab = "Time", ylab = "Dow Jones stock returns")

## test

set.seed(rseed)

djia_mt <- maxstat_test(djia_trafo ~ as.numeric(djia_time), distribution = approximate(9999), minprob = 0.1)

## maximally selected statistics

apply(abs(statistic(djia_mt, "standardized")), 2, max)

## critical values

qperm(djia_mt, 0.95)

## breakpoint

djia_point <- djia_time[floor(length(djia) * 0.1) + which.max(abs(statistic(djia_mt, "standardized")[,2]))]

## autocorrelation

djia_pre  <- window(djia, end = djia_point)

djia_post <- window(djia, start = djia_point + 1)

ar1 <- function(x, digits = 3) {

  x <- coredata(x)

  x <- x - mean(x)

  c(acf(x, plot = FALSE)$acf[2], acf(x^2, plot = FALSE)$acf[2])

}

ar1(djia_pre)

ar1(djia_post)

#######################

## Economic journals ##

#######################

## data

data("Journals", package = "AER")

Journals$age <- 2000 - Journals$foundingyear

Journals <- Journals[order(Journals$age),]

## model in root node

jour_lm  <- lm(log(subs) ~ log(price/citations), data = Journals)

## test in root node

set.seed(rseed)

(jour_mt <- maxstat_test(estfun(jour_lm) ~ age, data = Journals, distribution = approximate(9999), minprob = 0.1))

## test information

jour_critval <- qperm(jour_mt, 0.95^(1/4))

jour_process <- statistic(jour_mt, "standardized")

jour_process <- zoo(jour_process, as.numeric(sapply(strsplit(rownames(jour_process), "x <= "), tail, 1)))

colnames(jour_process) <- c("Intercept", "Slope")

jour_point <- time(jour_process)[apply(coredata(abs(jour_process)), 2, which.max)]

## visualization of test in root node

mypanel <- function(x, y, subscripts, groups, panel = panel.xyplot,

  col = 1, type = "p", pch = 20, lty = 1, lwd = 1, ...)

{

  col <- rep(as.list(col), length = nlevels(groups))

  type <- rep(as.list(type), length = nlevels(groups))

  pch <- rep(as.list(pch), length = nlevels(groups))

  lty <- rep(as.list(lty), length = nlevels(groups))

  lwd <- rep(as.list(lwd), length = nlevels(groups))

  for(g in 1:nlevels(groups)) {

    idx <- g == groups[subscripts]

    if (any(idx)) panel(x[idx], y[idx], ...,

      col = col[[g]], type = type[[g]], pch = pch[[g]],

      lty = lty[[g]], lwd = lwd[[g]])

    grid::grid.lines(y = grid::unit(0, "native"), gp = grid::gpar(col = "gray"))

    grid::grid.lines(y = grid::unit(jour_critval, "native"), gp = grid::gpar(col = 2))

  }

}

print(xyplot(abs(jour_process), panel = mypanel,

  xlab = "Time", type = "b", ylim = c(0, 6), as.table = TRUE))

### fit node with h(Y) = score and g(X) = maxstat_trafo

ytrf <- function(data) {

    ret <- estfun(lm(data[[1]] ~ data[[2]]))

    attr(ret, "assign") <- 1:2

    ret

}

xtrf <- function(data) trafo(data, numeric_trafo = maxstat_trafo)

mynode <- function(data) {

  vars <- c("society", "pages", "charpp", "age")

  sapply(vars, function(v) {

    f <- as.formula(paste("log(subs) + log(price/citations) ~", v))

    it <- independence_test(f, data = data,

      ytrafo = ytrf, xtrafo = xtrf, distribution = approximate(9999))

    c(statistic(it), 1 - (1 - pvalue(it))^length(vars))

  })

}

## fit all tree elements

jour_node <- factor(Journals$age <= 18, levels = c(TRUE, FALSE), labels = c("Node 2", "Node 3"))

jour_lm2 <- lm(log(subs) ~ log(price/citations), data = Journals[jour_node == "Node 2",])

jour_lm3 <- lm(log(subs) ~ log(price/citations), data = Journals[jour_node == "Node 3",])

## conduct tests in all leaves

set.seed(rseed)

jour_tree <- list(

  mynode(Journals),

  mynode(Journals[jour_node == "Node 2",]),

  mynode(data = Journals[jour_node == "Node 3",])

)

## fitted models

plot(log(subs) ~ log(price/citations), data = Journals,

  xlab = "log(price/citations)", ylab = "log(subscriptions)",

  pch = c(24, 21)[jour_node], bg = hcl(c(0, 240), 50, 70)[jour_node])

abline(coef(jour_lm2), col = hcl(  0, 80, 30), lty = 5, lwd = 1.7)

abline(coef(jour_lm3), col = hcl(240, 80, 30), lty = 1, lwd = 1.7)

legend("bottomleft", c(expression(age > 18), expression(age <= 18)),

  pch = c(19, 17), lty = c(1, 5), col = hcl(c(240, 0), 80, 30), bty = "n")