[R] faster GLS code

[Carlo Fezzi]

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