"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)
}