Revision d90516ef8fb8d7c848eb685f443359b2a267a100 authored by Roger Koenker on 18 September 2004, 00:00:00 UTC, committed by Gabor Csardi on 18 September 2004, 00:00:00 UTC
1 parent 926eb8c
khmal.R
"khmaladze.test" <-
function(fit, nullH = "location-scale", trim = c(0.25, 0.75) )
{
# Compute Khmaladze test statistics: fit is produced by rqProcess()
J <- standardize(fit, nullH)
Vtilde <- khmaladzize(fit$taus, fit$coef[1,], J$Vhat, nullH)
vtilde <- khmaladzize(fit$taus, fit$coef[1,], J$vhat, nullH)
trim <- ( fit$taus >= trim[1] & fit$taus <= trim[2] )
if(nullH == "location-scale") {
Tvtilde <- (vtilde-vtilde[,2])/sqrt(max(fit$taus)-min(fit$taus))
Tn <- max( apply(abs(Vtilde-Vtilde[,2])/
sqrt(max(fit$taus)-min(fit$taus)), 2,"sum")[trim] )
THn <- apply(abs(Tvtilde[,trim]),1,max)
}
else if( nullH == "location"){
Tn <- max( apply( abs(Vtilde-Vtilde[,2])/
sqrt(max(fit$taus)-min(fit$taus)), 2,"sum")[trim] )
THn <- apply(abs((vtilde[,trim]-vtilde[,2])/
sqrt(max(fit$taus)-min(fit$taus))),1,max)
}
x <- list(nullH = nullH, Tn = Tn, THn = THn)
class(x) <- "khmaladze"
x
}
"standardize" <-
function (fit, nullH = "location-scale")
{
Vhat <- fit$coefs
vhat <- fit$coefs
p <- nrow(vhat)
if (nullH == "location-scale") {
b <- matrix(0,2,p)
ER <- fit$coefs
for (i in 2:p) {
er <- lm(fit$coefs[i, ] ~ fit$coefs[1, ])
b[, i] <- er$coef
ER[i, ] <- er$resid
}
for (j in 1:length(fit$taus)) {
V <- fit$Hinv[, , j] %*% fit$J %*% fit$Hinv[, , j]
v <- V[-1, -1] + V[1, 1] * outer(b[2, -1], b[2, -1]) -
outer(V[-1, 1], b[2, -1]) - t(outer(V[-1, 1], b[2, -1]))
v <- solve(v)
v <- chol(0.5 * (v + t(v)))
Vhat[-1, j] <- v %*% ER[-1, j]
for (i in 2:dim(V)[1]) {
v <- V[i, i] + V[1, 1] * b[2, i]^2 - 2 * V[i, 1] * b[2, i]
vhat[i, j] <- ER[i, j]/sqrt(v)
}
}
}
else if (nullH == "location") {
b <- apply(fit$coefs, 1, mean)
ER <- fit$coefs - b
for (j in 1:length(fit$taus)) {
V <- fit$Hinv[, , j] %*% fit$J %*% fit$Hinv[, , j]
v <- V[-1, -1]
v <- solve(v)
v <- chol(0.5 * (v + t(v)))
Vhat[-1, j] <- v %*% ER[-1, j]
for (i in 2:dim(fit$Hinv)[1]) {
v <- V[i, i]
vhat[i, j] <- ER[i, j]/sqrt(v)
}
}
}
else stop(paste("unrecognized nullH: ",nullH))
list(Vhat=Vhat, vhat=vhat)
}
"khmaladzize" <-
function(tau, atau, Z, nullH)
{
p <- diff( tau )
p <- c( p[1], p )
score <- akj(atau, atau, p)
L <- length(tau)
gdot2 <- -score$psi
gdot <- cbind(rep(1,L),gdot2)
if(nullH == "location-scale")
{
gdot3 <- gdot2 * atau
gdot <- cbind(gdot, gdot3)
}
kmin <- 0
p <- nrow(Z)
# the v process
for (i in 1:p)
{
v <- Z[i,]
dtau <- diff(tau)
dtau <- c(dtau[1], dtau)
dv <- c(0,diff(v))
x1 <- gdot*sqrt(dtau)
x1 <- x1[L:1,]
y1 <- rev(dv/sqrt(dtau))
bhat <- lm.fit.recursive(x1,y1,int=FALSE)
bhat <- bhat[,L:1]
dvhat <- diag(gdot%*%bhat)*dtau
vhat <- cumsum(dvhat)
v <- v[kmin:L-kmin]
vhat <- vhat[kmin:L-kmin]
Z[i,] <- v-vhat
}
return(Z)
}
"rqProcess" <-
function(formula, data, taus=seq(0.2,0.8,by=0.1) )
{
z <- summary(lm(formula, data=data))
ols <- z$coef
vars <- names(z$coef[,1])
p <- length(z$coefficients[,1])
J <- solve(z$cov.unscaled)
#
# Compute RQ process
#
ntaus <- length( taus )
coefs<- matrix( 0, p, ntaus )
Hinv <- array( 0, c(p, p, ntaus) )
cat("taus: ")
for(i in 1:ntaus)
{
cat(taus[i]," ")
z <- summary(rq(formula,tau=taus[i],method="fn", data=data), se = "nid", hs=FALSE)
coefs[,i] <- z$coefficients[,1]
Hinv[,,i] <- z$Hinv
}
dimnames(coefs) <- list(vars,NULL)
x <- list(taus = taus, coefs=coefs, J=J, Hinv=Hinv)
class(x) <- "rqProcess"
x
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...