// This is a short bit of code to estimate a linear regression model that allows // the residual variance (actually the natural log of the standard deviation) to // be a function of a set of predictor variables. The model is estimated by // maximum likelihood. // // This sort of model is discussed pp 28-34 of Scott. R. Eliason (1993) // Maximum Likelihood Estimation: Logic and Practice, Sage. // Eliason calls it the heteroscedastic normal pdf model. // // Bruce Western and Deirdre Bloome discuss various ways of doing this sort of // thing in their 2009 Sociological Methodology article 'Variance Function // Regressions for Studying Inequality' // // Western and Bloome give the Stata code for one of their methods at the end of // their article. You should note that there is a trivial typo in the printed // version of their code that will generate an error message. I've corrected // the typo and give their code below mine. I find my way of doing things much // more direct, but I haven't investigated how it performs in very many cases // and there may be advantages to their method which I haven't appreciated. In // any case if you use the Western/Bloome code you should acknowledge them and // not me! // Begin by generating some data. A binary variable x1, a continuous variable x2 // and a response variable y1. The important feature of the data is that // the residual variation is much greater if x1 = 1 than if x1 = 0. This is // what we want to recover from the model. set obs 10000 gen z = uniform() gen x1 = 0 if z <0.5 replace x1 =1 if z >= 0.5 gen x2 = rnormal(0,5) gen y1 = 1.5*x1 + 0.5*x2 + rnormal(0,1) if x1==0 replace y1 = 1.5*x1 + 0.5*x2 + rnormal(0,5) if x1==1 // Now establish the target tab x1, summ(y1) // Note what the standard deviation is if x1 = 0 and if x1 = 1 // Now define the likelihood to be maximised and the arguments capture program drop hetnormal program define hetnormal version 14.1 args lnf mu lnsigma quietly replace `lnf' = /// ln(normalden( \$ML_y1, `mu', exp(`lnsigma'))) end // This bit defines the model and tells stata to maximise the likelihood // We'll start with just one predictor variable - x1 ml model lf hetnormal /// (mu: y1 = x1 ) /// (lnsigma: x1), /// diparm (lnsigma, exp label("sigma")) ml check ml max // The slopes for the means are reported in the panel labeled mu. // The slopes for the log standard deviation are reported in the panel lnsigma // // sigma is exp(_cons) ie the standard deviation if x1 = 0 // The standard deviation if x1 = 1 is exp(x1 + _cons) // This is the Western/Broome code // NB the scale for the model of the residual variation is log of the variance // ie to compare the results with those produced by my code you need to // exponentiate and take the square root reg y1 x1 predict R, r gen R2=R^2 glm R2 x1, family(gamma) link(log) predict S2, mu gen LOGLIK = (-.5*(ln(S2)+(R2/S2) )) egen LL0 = sum(LOGLIK) di LL0 * Updating beta and lambda coefficients; gen DLL=1 while DLL > .00000001 { drop R quietly: reg y1 x1 [aw=1/S2] drop S2 predict R, r replace R2=R^2 est store BETA quietly: glm R2 x1, family(gamma) link(log) predict S2, mu est store LAMBDA replace LOGLIK = (-.5*(ln(S2)+(R2/S2))) egen LLN = sum(LOGLIK) di LLN replace DLL=LLN-LL0 replace LL0=LLN drop LLN } est table BETA LAMBDA, b se // This adds another predictor which is continuous x2. It should have no effect on residual variation ml model lf hetnormal /// (mu: y1 = x1 x2) /// (lnsigma: x1 x2), /// diparm (lnsigma, exp label("sigma")) ml check ml max