Variance Function Regressionをmlexpで解く

Western and Bruceが示しているコードではなく，尤度関数を最大化する方法をとる．
これはstataのmlexpを用いた方法である．

このブログを参考にする．
ただし，ローカルマクロに関してはwindowsで利用可能なように書き換えた．

sysuse auto, clear
local x = "weight mpg foreign"
local z = "weight mpg foreign"
reg price `x'
predict R, r
gen R2=R^2

glm R2 `z', family(gamma) link(log)
predict S2, mu

gen LOGLIK = -(1/2)*(ln(S2)+(R2/S2))
egen LL0 = sum(LOGLIK)
display LL0
gen DLL=1
while DLL > .00001 {
drop R
quietly reg price `x' [aw=1/S2] 
drop S2
predict R, r
replace R2=R^2
est store BETA
quietly glm R2 `z', family(gamma) link(log)
predict S2, mu
est store LAMBDA
replace LOGLIK = -(1/2)*(ln(S2)+(R2/S2))
egen LLN = sum(LOGLIK)
display LLN
replace DLL=LLN-LL0 
replace LL0=LLN
drop LLN
}
est tab BETA LAMBDA, b se stat(N r2)

上記のモデルをmlexpで解くときは以下のようになる．
尤度関数については，ここの13〜14ページを参照．

mlexp(-0.5*(({l1}*weight+{l2}*mpg+{l3}*foreign+{l0})+(price-({b1}*weight+{b2}*mpg+{b3}*foreign+{b0}))^2*exp(-({l1}*weight+{l2}*mpg+{l3}*foreign+{l0}))))

少し異なる結果が得られた．mlexpで解くと標準誤差が少し大きく推定される．
もともとのコードでは以下のような結果が得られる．

----------------------------------------
    Variable |    BETA        LAMBDA    
-------------+--------------------------
_            |
      weight |   .7621136               
             |  .39582551               
         mpg | -7.3244202               
             |   40.84604               
     foreign |  1022.2735               
             |  335.75554               
       _cons |  2821.8683               
             |  1906.8068               
-------------+--------------------------
R2           |
      weight |               .00220443  
             |               .00047787  
         mpg |               .00479197  
             |               .06017879  
     foreign |               2.0162373  
             |               .48658851  
       _cons |               7.4684108  
             |               2.6853067  
-------------+--------------------------
Statistics   |                          
           N |         74           74  
          r2 |  .20506016               
----------------------------------------
                            legend: b/se

mlexpを用いると以下の結果が得られる．

Maximum likelihood estimation

Log likelihood = -585.56457                     Number of obs     =         74

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         /l1 |   .0022048   .0004572     4.82   0.000     .0013088    .0031009
         /l2 |   .0047836   .0549507     0.09   0.931    -.1029178    .1124851
         /l3 |   2.016866   .4955487     4.07   0.000     1.045609    2.988124
         /l0 |   7.467199   2.482776     3.01   0.003     2.601048    12.33335
         /b1 |   .7614853   .4257266     1.79   0.074    -.0729235    1.595894
         /b2 |   -7.30936   41.04674    -0.18   0.859    -87.75949    73.14077
         /b3 |   1021.795   364.0346     2.81   0.005     308.3001     1735.29
         /b0 |   2822.992   1927.694     1.46   0.143    -955.2195    6601.203
------------------------------------------------------------------------------

対数尤度も少し違うのは何に起因するのだろうか．