[R] faster GLS code

[Carlo Fezzi]

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: rvaradhan at jhmi.edu
>>
>>
>> ----- Original Message -----
>> From: Carlo Fezzi <c.fezzi at uea.ac.uk>
>> Date: Thursday, January 7, 2010 12:13 pm
>> Subject: [R] faster GLS code
>> To: r-help at 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.fezzi at 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 at r-project.org mailing list
>>>
>>>  PLEASE do read the posting guide
>>>  and provide commented, minimal, self-contained, reproducible code.
>>
>> ______________________________________________
>> R-help at 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:cberry at tajo.ucsd.edu	            UC San Diego
> http://famprevmed.ucsd.edu/faculty/cberry/  La Jolla, San Diego 92093-0901
>
>
>