# How to implement a variance function in R?

I want to estimate a multivariate variance function in R. That is, I want to allow the variance (as well as the mean) to vary according to some set of independent variables.

In this particular case, I want to estimate the effects of a set of typical demographic covariates (age, race, education) on the variance of logged wages.

What is a good way to implement this in R? Is there a package that simplifies this?

It may be that this is only a search away - but having searched on the R help pages, Google, Rseek, and StackOverflow, I can't find anything relevant under "variance function" or similar.

Thanks for your responses -- I will try to clarify my question.

I am working in a maximum likelihood framework. I can code this by hand from the log-likelihood, but the real data set has a lot of variables and "optim" is very slow, so I would like to find a package in R that makes this more computationally efficient.

I start with the log-likelihood for a basic OLS regression: $$\text{ln }L = \sum (-\frac{1}{2} (\text{ln }\sigma^2 - \frac{(y - xB)^2}{\sigma^2}))$$ Then I relax the assumption of constant variance (homoskedasticity) and redefine the variance as: $$\sigma^2 = exp(Z*\gamma)$$ where $Z$ is the matrix of variables affecting $\sigma^2$. (Exponentiate so that you don't end up with $\sigma^2$ less than zero.) When I substitute the reparameterization of $\sigma^2$ into the original log-likelihood and code the new log-likelihood function in R, I get this:

ll.normal.vary <- function (par, X, Y, Z) {
beta  <-par[1:ncol(X)]
gamma <- par[(ncol(X)+1):(ncol(X)+ncol(Z))]
-1/2* sum((Z %*% gamma) + ((Y - X %*% beta)^2)/exp(Z %*% gamma))
}


Then I optimize:

v.optim1 <- optim (par = start1, fn=ll.normal.vary, X=x.mat, Y=y.vec, Z=z.mat,
method = "BFGS", hessian = F, control = list(fnscale = -1))
v.optim1$par v.optim1$value


Here are some sample data if you want to test it:

var1   <- c(0,0,0,1,1,0)
var2   <- c(.28, .07, -.05, .38, .08, -.1)
var3   <- c(-.11, -.17, -.17, -.05, .1, -.01)
x.mat  <- cbind(var1, var2, var3)
y.vec  <- c(.46, .77, .49, .59, .60, .44)
z.mat  <- cbind(var1, var2)
start1 <- rep(0.1, ncol(x.mat)+ncol(z.mat))


Thanks again for any tips.

-
There are lots of statistical frameworks for estimating variances, ranging from maximum likelihood estimation to complex dynamical/time-series approaches like GARCH ... Once you know what statistical framework you want to use, and have read a little bit to understand how it works, you can figure out what R's capabilities are (e.g. asking here again, or using library(sos); findFn("whatever") once you know what keywords you are looking for) –  Ben Bolker Apr 4 '12 at 21:15
I'm not sure, your second paragraph let me think of "variance components", which can be done in R with the lmer function from the lme4 package. –  psj Apr 5 '12 at 14:40
@boston192 I presume you are the author of the question; if so the system failed to merge you. In this case, please go to SO, login there (as user1257313) and come back here. This should correct the ownership problem. –  mbq Apr 10 '12 at 14:26
@mbq Thanks for the response, I've re-logged in under my SO login. Hope that sorts it out. –  user1257313 Apr 10 '12 at 14:31
@user1257313 Great, thanks! –  mbq Apr 10 '12 at 14:33

## migrated from stackoverflow.comApr 5 '12 at 12:30

This question came from our site for professional and enthusiast programmers.