time format in r axis

x=c(paste(“81m”, 1:12, sep=”"),paste(“82m”, 1:12, sep=”"))
ind=1:24
y=cumsum(rnorm(24,0))
par(lwd=2)
plot(y~ind,type=’b',xaxt=’n',col=4,ylab=”,xlab=”);grid()

axis(1,seq(1,24,2),labels=F,tcl=-0.3)
#down axis
axis(1,seq(2,24,2),labels=F,tcl=-0.5,lwd.ticks=2)
mtext(x[seq(1,24,2)],1,at=seq(1,24,2),padj=3)
#down axis names
mtext(x[seq(2,24,2)],1,at=seq(2,24,2),padj=1)

Speed test

%Matlab

m=ones(200000,10);
ret=0;
tic
for i=1:20,
ret=ret+sum(sum(m,1),2);
end
toc

/*Gauss*/
m=ones(200000,10);
ret=0;

starttime=date;
for i (1,20,1);
ret=ret+sumc(sumr(m));
endfor;
endtime=date;
print etstr( ethsec(starttime,endtime) );

/*OX*/
#include
main(){
decl tim,i,m,ret,k;
m = ones(200000,10);
ret =0;

tim = timer();
for (i=0;i<20;++i){
ret=ret+sumc(sumr(m));
}
print(“time lapsed: “, timespan(tim, timer()), “\n”);
}

#R
# File name: loop_rowSums.r

my.loop <- 20
m.dim <- list(nrow = 200000, ncol = 10)
m <- matrix(1, nrow = m.dim$nrow, ncol = m.dim$ncol)
ret <- 0

start <- Sys.time()
for(k in 1 : my.loop){
ret <- ret + sum(rowSums(m))
}
Sys.time() – start
==========================

Speed:

Matlab(0.119)
OX(0.2)
Gauss(0.14)
R(1.625)

http://www.math.ncu.edu.tw/~chenwc/R_note/index.php?item=loop

R barplot

 x1 <- c(23.2,34.5,76.3,65.8,12.6)
  x2 <- c(15.6,12.4,21.8,20,5.2)
  A <- gl(5,1,5,labels=c("a1","a2","a3","a4","a5"))
  data <- cbind(x1,x2)
  rownames(data) <- levels(A)
  barplot(x1,names.arg=levels(A))
  barplot(t(data),beside=T,ylim=c(0,100),legend.text=colnames(data),
    col=c("grey50","grey80"),ylab="Fréquence")

Simulate price process

http://www.kevinsheppard.com/matlab/Assignment2.html
use R and OX

probit ml in stata and R

/* Stata */

clear
set obs 1000
gen x=invnorm(uniform())
gen e=invnorm(uniform())
gen y=0.5+0.4*x+e
replace y=0 if y<=0
replace y=1 if y>0

ols and probit diff

n<-10000
disper<-0.6
x<-rnorm(n)
e<-rnorm(n,disper)
a<-0.3;b<-0.6;
y<-a+b*x+e
y<-ifelse(y<0,0,1)
g<-y~x
summary(lm(g))
summary(glm(g,family=binomial(link=probit)))

Title:
Results on the bias and inconsistency of ordinary least squares for the linear probability model

Author
William C. Horrace^a, and Ronald L. Oaxaca^b

Limitations of the Linear Probability Model (LPM) are well-known. OLS estimated probabilities are not bounded on the unit interval, and OLS estimation implies that heteroscedasticity exists. Conventional advice points to probit or logit as the standard remedy, which bound the maximum likelihood estimated probabilities on the unit interval. However, the fact that consistent estimation of the LPM may be difficult does not imply that either probit or logit is the correct specification of the probability model; it may be reasonable to assume that probabilities are generated from bounded linear decision rules. Theoretical rationalizations for the LPM are in Rosenthal (1989) and Heckman and Snyder, 1977 J.J. Heckman and J.M. Snyder Jr., Linear probability models of the demand for attributes with an empirical application to estimating the preferences of legislators, Rand Journal of Economics 28 (1977), pp. S142–S189.Heckman and Snyder (1977).

Despite the attractiveness of logit and probit for estimating binary dependent variable models, OLS on the LPM is still used. Recent applications include Klaassen and Magnus (2001), Bettis and Fairlie (2001), Lukashin (2000), McGarry (2000), Fairlie and Sundstrom (1999), Reiley (2005), and Currie and Gruber (1996). Empirical rationales for the LPM specification are plentiful. McGarry appeals to ease of interpretation of estimated marginal effects, while Reiley cites a perfect correlation problem associated with the probit model. Fairlie and Sundstrom prefer LPM because it implies a simple expression for the change in unemployment rate between two censuses. Bettis and Farlie choose LPM because of an extremely large sample size and other simplifications implied by it. Lukashin uses the LPM, because it lends itself to a model selection algorithm based on an adaptive gradient criterion. Currie and Gruber state that logit, probit, and OLS are similar for their data and only report LPM results.

Other rationales for the OLS on the LPM are complications of probit/logit models in certain contexts. Klaassen and Magnus cite panel data complications in their tennis example and select OLS. OLS is perhaps justified in simultaneous equations/instrumental variable methods. The presence of dummy endogenous regressors is problematic if the DGP is assumed to be probit or logit; these problems were first considered by Heckman (1978 ). While perhaps less popular than logit and probit, OLS on the LPM model still finds its way into the literature for various reasons.

Some well-known LPM theorems are provided in Amemiya (1977). Econometrics textbooks (e.g., Greene, 2000), acknowledge complications leading to biased and inconsistent OLS estimates. Nevertheless, the literature is not clear on the precise conditions when OLS is problematic. This note rigorously lays out these conditions, derives the finite-sample and asymptotic biases of OLS, and provides additional results that highlight the appropriateness or inappropriateness of OLS estimation of the LPM. Finally, we suggest a trimmed sample estimator that could reduce OLS bias.

svg in R

library(“RSVGTipsDevice”)
library(“SemiPar”)
data(swiss)
setwd(‘c:/data/’)
devSVGTips(“swiss.svg”, height=7, width=10, toolTipMode=1,
title=”Swiss  catholic  and agriculture  plot”)
par(font.lab=8,font.main=9,cex=1.5)
plot(swiss$Agriculture~swiss$Catholic,ylim=c(0,130),
type=”n”, xlab=”catholic”,ylab=”Agri”, main=”catholic and agriculture”)
grid()
cnames<-rownames(swiss)
swiss.1<-cbind(cnames,swiss)
for (i in seq(len=nrow(swiss.1))) {
setSVGShapeToolTip(title=swiss.1[i,"cnames"],
desc=paste(” agri=”, swiss.1[i,"Agriculture"],” ,catholic”, swiss.1[i,"Catholic"]))
points(swiss.1[i,"Agriculture"], swiss.1[i,"Catholic"], pch=12,
cex=2, col=ifelse(swiss.1[i,"Catholic"]>50,”blue”,”red”))
}
dev.off()

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”)