robustreg in R(not complete)

packages hccm

robreg<-function(y,x){
n<-length(y)
k<-dim(cbind(1,x))[2]
xpx<-crossprod(x)
ypy<-crossprod(x,y)
beta<- solve(xpx)%*%ypy
#b<- solve(qr(x),y)
#b<- solve(svd(x),y)
e <- y-x%*%b
s2 <- crossprod(e)/(n-k)
vbhom <- s2[1]*solve(xpx)
sebhom <- sqrt(diag(vbhom))
tbhom <- b/sebhom
pbhom <- 2*(1-pt(abs(tbhom),n-k))
# OLS standard errors for heteroskedastic errors
esq <- e^2
Vbhet <- solve(xpx)%*%(t(x)%*%(x*esq))%*%solve(xpx)%*%(n/(n-k))
sebhet <- sqrt(diagonal(Vbhet))
tbhet <- b/sebhet
pbhet <- 2*(1-pt(abs(tbhet),n-k))
}

summaryw <- function( model) {
s <- summary( model)
X <- model.matrix( model)
u2 <- residuals( model)^2
XDX <- 0
## here one needs essentially to calculate X’DX.  But due to the fact that D
## is huge (NxN), it is better to do it with a cycle.
for( i in 1:nrow( X)) {
XDX <- XDX + u2[i]*X[i,]%*%t(X[i,])
}
XX1 <- solve( t(X)%*%X)
varcovar <- XX1 %*% XDX %*% XX1
stdh <- sqrt( diag( varcovar))
t <- model$coefficients/stdh
p <- 2*pnorm( -abs( t))
results <- cbind( model$coefficients, stdh, t, p)
dimnames(results) <- dimnames( s$coefficients)
results
}

mata
// Create y vector and X matrix from Stat data set using st_view command
st_view(y=., ., “wage76″)
st_view(X=., ., (“grade76″, “age76″, “cons”))
Xnames = (“grade76″, “age76″, “cons”)   // used to label output
// Compute OLS estimator
b = invsym(X’X)*X’y
// Compute OLS standard errors for homoskedastic errors
e = y – X*b
n = rows(X)
k = cols(X)
s2 = (e’e)/(n-k)
Vbhom = s2*invsym(X’X)
sebhom = sqrt(diagonal(Vbhom))
tbhom = b:/sebhom
pbhom = 2*ttail(n-k, abs(tbhom))
// OLS standard errors for heteroskedastic errors
esq = e:*e
Vbhet = invsym(X’X)*(X’(X:*esq))*invsym(X’X)*(n/(n-k))
sebhet = sqrt(diagonal(Vbhet))
tbhet = b:/sebhet
pbhet = 2*ttail(n-k, abs(tbhet))
// Print results
(b, sebhom, tbhom, pbhom)
(b, sebhet, tbhet, pbhet)
// Easier to read output includes regressor names
(Xnames’, strofreal(b), strofreal(sebhom), strofreal(tbhom), strofreal(pbhom))
(Xnames’, strofreal(b), strofreal(sebhet), strofreal(tbhet), strofreal(pbhet))
end

* CHECK: OLS using regress
reg wage76 grade76 age76
reg wage76 grade76 age76, robust