R computing speed

Dear helpers,

I am using R version 2.5.1 to estimate a multinomial logit model using my
own maximum likelihood function (I work with share data and the default
function of R cannot deal with that).

However, the computer (I have an Athlon XP 3200+ with 512 GB ram) takes
quite a while to estimate the model.

With 3 categories, 5 explanatory variables and roughly 5000 observations it
takes 2-3 min. For 10 categories and 10 explanatory variables (still 5000
obs) more than 1 hour.

Is there any way I can speed up this process? (Modifying the code or
modifying some R options maybe?)

I would be really grateful if anybody could help me with this issue, I
attach my code below.

Many thanks,

Carlo

***************************************
Carlo Fezzi

Centre for Social and Economic Research
on the Global Environment (CSERGE),
School of Environmental Sciences,
University of East Anglia,
Norwich, NR4 7TJ
United Kingdom.

***************************************



# MULTILOGIT

# This function computes the estimates of a multinomial logit model

# inputs:       a matrix vector of 1 and 0 (y) or of shares
#               a matrix of regressors (x) - MUST HAVE COLUMN NAMES! -
#               names of the variables, default = colnames(x)
#               optimization methods, default = 'BFGS'
#               base category, default = 1
#               restrictions, default = NULL
#               weights, default all equal to 1


# outputs:      an object of class "multilogit.c"

# McFadden D. (1974) "Conditional logit analysis of qualitative choice
behavior", in Zarembka P. (ed.), Frontiers in Econometrics, Academic Press.


multilogit.c <- function(y, xi, xi.names = colnames(xi), c.base=1,
rest=NULL, w = rep(1,nrow(y)), method='BFGS')
{

       n.obs <- sum(w)
       xi<-cbind(1,xi)
       colnames(xi)[1]<-"Intercept"

       nx<-ncol(xi)
       ny<-ncol(y)

       beta<-numeric(nx*ny)

       negll<- function(beta,y,xi)
       {
               beta[rest]<-0
               beta[(((c.base-1)*nx)+1):(c.base*nx)]<-0
               lli <- y  * (xi%*%matrix(beta,nx,ny) - log ( apply(exp(
xi%*%matrix(beta,nx,ny)) ,1,sum ) )     )
               lli<-lli*w
               -sum(lli)
       }

       pi<- apply((y*w),2,mean)/mean(w)

       ll0 <- (t(pi)%*%log(pi))*sum(w)

       result<-c(      optim(par = rep(0,nx*(ny)), fn = negll, y=y, xi=xi,
hessian=T, method=method),
                       list(varnames=xi.names, rest=rest, nx=nx, ny=ny,
npar=nx*(ny-1)-length(rest), ll0=ll0,   pi=pi, xi=xi,
n.obs=n.obs,c.base=c.base,w=w))

       result$par <- result$par[-(((c.base-1)*nx)+1):-(c.base*nx)]
       result$hessian <-
result$hessian[-(((c.base-1)*nx)+1):-(c.base*nx),-(((c.base-1)*nx)+1):-(c.ba
se*nx)]

       class(result)<-"multilogit.c"
       return(result)
}