"htiv" <- function(x, y, id, d, z = NULL) { #input: # x design matrix partitioned as given in d # y response vector # id strata indicator # d list of column numbers indicating partitioning of x # x[,d[[1]]] is x1 -- exogonous time varying vars # x[,d[[2]]] is x2 -- endogonous time varying vars # x[,d[[3]]] is z1 -- exogonous time invariant vars # x[,d[[4]]] is z2 -- endogonous time invariant vars # z may contain excluded exogonous variables if there are any # NB. intercept is automatically included x <- as.matrix(cbind(x, 1)) Tx <- PQ(x, id) d[[3]] <- c(d[[3]], dim(x)[2]) Z <- cbind(z, Tx\$Ph[, d[[1]]], Tx\$Qh[, d[[1]]], Tx\$Qh[, d[[2]]], x[, d[[ 3]]]) r <- tsls(x, Z, y, int = F)\$resid Ti <- table(id) Ti.inv <- 1/table(id) rdot2 <- tapply(r, id, mean)^2 v <- lm(rdot2~Ti.inv) v <- v\$coef theta <- as.vector(sqrt(v[2]/(v[2] + v[1] * Ti[Tx\$is]))) x <- x - (1 - theta) * Tx\$Ph y <- y - (1 - theta) * PQ(y, id)\$Ph fit <- tsls(x, Z, y, int = F) list(fit=fit,v=v) } ______________________________Cut Here ________________________________ "PQ" <- function(h, id) { if(is.vector(h)) h <- matrix(h, ncol = 1) Ph <- unique(id) Ph <- cbind(Ph, table(id)) for(i in 1:ncol(h)) Ph <- cbind(Ph, tapply(h[, i], id, mean)) is <- tapply(id, id) Ph <- Ph[is, - (1:2)] Qh <- h - Ph list(Ph=as.matrix(Ph), Qh=as.matrix(Qh), is=is) } ______________________________Cut Here ________________________________ "tsls" <- function(x, z, y, int = T) { # two stage least squares if(int){ x <- cbind(1, x) z <- cbind(1, z) } xhat <- lm(x~z-1)\$fitted.values R <- lm(y ~ xhat -1) R\$residuals <- c(y - x %*% R\$coef) return(R) }