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 ------------------------------------------------------------------------------
対数尤度も少し違うのは何に起因するのだろうか.