Dear Ravi and Chuck, thanks a lot for your suggestions. I should have make it clear earlier, but unfortunately I need "betav" because this is part of a longer code to estimate a bayesian SUR tobit. Ravi, in your approach is there are way to compute "betav" as well? Also, I was aware of the model reduction to a simple OLS when the regressors are the same, but I still need the full GLS estimator because of its use with the rest of the model.
Thanks a lot, Carlo > On Thu, 7 Jan 2010, Ravi Varadhan wrote: > >> Try this: >> >> >> X <- kronecker(diag(1,3),x) >> Y <- c(y) # stack the y in a vector >> >> # residual covariance matrix for each observation >> covar <- kronecker(sigma,diag(1,N)) >> >> csig <- chol2inv(covar) >> betam2 <- ginv(csig %*% X) %*% csig %*% Y >> >> This is more than 2 times faster than your code (however, it doesn't >> compute `betav') . >> > > Faster still (by a wide margin) if X will truly be of that form: > >> B <- coef(lm(y~0+.,as.data.frame(x))) >> all.equal( as.vector((B)), as.vector(betam)) > [1] TRUE > > When X is of that form, the covariance matrix drops out of the > computation. > > :) > > Chuck > > >> Here is a timing comparison: >> >> # Your method >> # GLS betas covariance matrix >> system.time({ >> inv.sigma <- solve(covar) >> betav <- solve(t(X)%*%inv.sigma%*%X) >> >> # GLS mean parameter estimates >> betam <- betav%*%t(X)%*%inv.sigma%*%Y >> }) >> >> # New method >> system.time({ >> csig <- chol2inv(covar) >> betam2 <- ginv(csig %*% X) %*% csig %*% Y >> }) >> >> all.equal(betam, betam2) >> >>> # GLS betas covariance matrix >>> system.time({ >> + inv.sigma <- solve(covar) >> + betav <- solve(t(X)%*%inv.sigma%*%X) >> + >> + # GLS mean parameter estimates >> + betam <- betav%*%t(X)%*%inv.sigma%*%Y >> + }) >> user system elapsed >> 1.14 0.51 1.76 >>> >>> system.time({ >> + csig <- chol2inv(covar) >> + betam2 <- ginv(csig %*% X) %*% csig %*% Y >> + }) >> user system elapsed >> 0.47 0.08 0.61 >>> >>> all.equal(betam, betam2) >> [1] TRUE >>> >> >> >> Hope this helps, >> Ravi. >> ____________________________________________________________________ >> >> Ravi Varadhan, Ph.D. >> Assistant Professor, >> Division of Geriatric Medicine and Gerontology >> School of Medicine >> Johns Hopkins University >> >> Ph. (410) 502-2619 >> email: rvarad...@jhmi.edu >> >> >> ----- Original Message ----- >> From: Carlo Fezzi <c.fe...@uea.ac.uk> >> Date: Thursday, January 7, 2010 12:13 pm >> Subject: [R] faster GLS code >> To: r-help@r-project.org >> >> >>> Dear helpers, >>> >>> I wrote a code which estimates a multi-equation model with generalized >>> least squares (GLS). I can use GLS because I know the covariance >>> matrix of >>> the residuals a priori. However, it is a bit slow and I wonder if >>> anybody >>> would be able to point out a way to make it faster (it is part of a >>> bigger >>> code and needs to run several times). >>> >>> Any suggestion would be greatly appreciated. >>> >>> Carlo >>> >>> >>> *************************************** >>> Carlo Fezzi >>> Senior Research Associate >>> >>> Centre for Social and Economic Research >>> on the Global Environment (CSERGE), >>> School of Environmental Sciences, >>> University of East Anglia, >>> Norwich, NR4 7TJ >>> United Kingdom. >>> email: c.fe...@uea.ac.uk >>> *************************************** >>> >>> Here is an example with 3 equations and 2 exogenous variables: >>> >>> ----- start code ------ >>> >>> >>> N <- 1000 # number of observations >>> library(MASS) >>> >>> ## parameters ## >>> >>> # eq. 1 >>> b10 <- 7; b11 <- 2; b12 <- -1 >>> >>> # eq. 2 >>> b20 <- 5; b21 <- -2; b22 <- 1 >>> >>> # eq.3 >>> b30 <- 1; b31 <- 5; b32 <- 2 >>> >>> # exogenous variables >>> >>> x1 <- runif(min=-10,max=10,N) >>> x2 <- runif(min=-5,max=5,N) >>> >>> # residual covariance matrix >>> sigma <- matrix(c(2,1,0.7,1,1.5,0.5,0.7,0.5,2),3,3) >>> >>> # residuals >>> r <- mvrnorm(N,mu=rep(0,3), Sigma=sigma) >>> >>> # endogenous variables >>> >>> y1 <- b10 + b11 * x1 + b12*x2 + r[,1] >>> y2 <- b20 + b21 * x1 + b22*x2 + r[,2] >>> y3 <- b30 + b31 * x1 + b32*x2 + r[,3] >>> >>> y <- cbind(y1,y2,y3) # matrix of endogenous >>> x <- cbind(1,x1, x2) # matrix of exogenous >>> >>> >>> #### MODEL ESTIMATION ### >>> >>> # build the big X matrix needed for GLS estimation: >>> >>> X <- kronecker(diag(1,3),x) >>> Y <- c(y) # stack the y in a vector >>> >>> # residual covariance matrix for each observation >>> covar <- kronecker(sigma,diag(1,N)) >>> >>> # GLS betas covariance matrix >>> inv.sigma <- solve(covar) >>> betav <- solve(t(X)%*%inv.sigma%*%X) >>> >>> # GLS mean parameter estimates >>> betam <- betav%*%t(X)%*%inv.sigma%*%Y >>> >>> ----- end of code ---- >>> >>> ______________________________________________ >>> R-help@r-project.org mailing list >>> >>> PLEASE do read the posting guide >>> and provide commented, minimal, self-contained, reproducible code. >> >> ______________________________________________ >> R-help@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-help >> PLEASE do read the posting guide >> http://www.R-project.org/posting-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> > > Charles C. Berry (858) 534-2098 > Dept of Family/Preventive > Medicine > E mailto:cbe...@tajo.ucsd.edu UC San Diego > http://famprevmed.ucsd.edu/faculty/cberry/ La Jolla, San Diego 92093-0901 > > > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.