mata poisson

From Colin Cameron
http://cameron.econ.ucdavis.edu/stata/stata.html

* This program does

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

optimize in R et stata

stata
version 10
mata
void myeval(todo, x, y, g, H){
y = exp(-x^2 + x – 3)
}

S = optimize_init()
optimize_init_evaluator(S, &myeval())
optimize_init_params(S, 0)
x = optimize(S)
x
end

R
fr<-function(x){
-exp(-x^2+x-3)}
optim(c(0),fr,method=”BFGS”)

R:two axis label

http://www.apsnet.org/education/AdvancedPlantPath/Topics/RModules/Doc1/01_Nematode.html
## Set up the vector of nematode populations
npop <- c(220,180,150,250,270,300,500,580,580,1000,380,100)

## Set up vector of associated months
names(npop) <-c (
        "Jan",
        "Feb",
        "Mar",
        "Apr",
        "May",
        "Jun",
        "Jul",
        "Aug",
        "Sep",
        "Oct",
        "Nov",
        "Dec"
)

## Set up the bar plots of the nematode populations,
try help(barplot) for more information
barplot(
        npop,
        axes=FALSE,
        main="Nematode  Population",
        ylim=c(0,1200),
        col="orange",
        xlab="Month",
        ylab="Nematodes per 100 cc soil")
        axis(2,seq(0,1200,by=200)
)

## Set up the new plot to be drawn on top of the barplot,
try help(par) for more information
par(new=T)

# This sets the average monthly temperature
# for the line above the population
y <-c (15,17,20,24,26,28,30,31,30,28,19,16)

# This positions the points on the X axis
# correctly in the middle over the top of the bar chart.
y2 <-c (0.5,1.45,2.45,3.5,4.5,5.35,6.5,7.6,8.5,9.55,10.65,11.5)

## Plot the temperature points
plot(
        y2,
        y,
        ylim=c(0,35),
        xlab=NA,
        ylab=NA,
        xlim=c(0,12),
        axes=FALSE,
        type="o",
        pch=22,
        col="mediumblue"
)

axis(
        4,
        seq(0,35,by=5),
        labels=FALSE
)

## Lay out the points for temperature
points(
        y2,
        y,
        pch=22,
        col="mediumblue"
)

## Label the right axis, try help(mtext) for more information
mtext(
        c("0C","5C","10C","15C","20C","25C","30C","35C"),
        4,
        at=seq(0,35,by=5),
        line=.25
)

## Label the right axis, try search(mtext) for more information
mtext(
        "Average Soil Temperature",
        4,
        line=1,
        cex=1
)

## Create a legend
legend(
        "topleft",
        c("Soil Temperature","Nematode Population"),
        pch=c(22,NA),
        lty=c(1,NA),
        col=c("mediumblue",NA),
        inset=0.01
)

## Create an orange rectangle next to
the Nematode Population label for the legend
rect(-0.2,32.4,0.4,33.4,col="orange")