faster GLS code

____________________________________________________________________

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.