Have you looked into the residual options in xtmixed?
residual_options are by(varname) and t(varname).
by(varname) is for use within the residuals() option and
specifies that a set of distinct residual-error parameters be
estimated for each level of varname. In other words, you use
by() to model heteroskedasticity.