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