Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: re: maximum likelihood problem


From   "Ben Spong" <bspong.stata@googlemail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: re: maximum likelihood problem
Date   Sun, 9 Dec 2007 11:43:22 +0000

Thanks for the help Kit! It allowed me to change a bit of the program
and almost have it correct. The program now seems to is the following:

program mymlprog
	version 10.0
	args lnf mu1 mu2 sigma

	quietly replace `lnf' = log(1/sqrt(2*_pi*`sigma'^2))-
(1/(2*`sigma'^2))*($ML_y1+`mu1')^2 if $ML_y1 <0

	quietly replace `lnf' =log(normal((`mu2')/`sigma')-
normal((`mu1')/`sigma')) if $ML_y1 == 0

	quietly replace `lnf' =log(1/sqrt(2*_pi*`sigma'^2))-
(1/(2*`sigma'^2))*($ML_y1+`mu2')^2 if $ML_y1 > 0

end

ml model lf mymlprog (mu1: ret=mktret) (mu2: ret=mktret) /sigma
ml check
ml maximize

Although the results are very similar to the ones I get in Excel, they
aren't quite the same because the beta coefficients different in the
two equations. However, I need them to be the same... Is there a way
of forcing Stata to produce these results?

Best,

B.

On Dec 8, 2007 4:19 PM, Kit Baum <baum@bc.edu> wrote:
> Ben wrote
>
> The program I've written in Stata is the following
> (which, for obvious reasons, isn't working):
>
> program mymlprog
>         version 10.0
>         args y x
>
>         quietly replace `y' = ln(1/sqrt(2*_pi*Sigma_hat^2))-
> ((1/(2*Sigma_hat^2))*(`y'+Alpha_1-Beta_hat*`x')^2 if $ML_y1 <0
>
>         quietly replace `y' =
> ln(normal((Alpha_2-Beta_hat*`x')/Sigma_hat)-normal((Alpha_1-
> Beta_hat*`x')/Sigma_hat))
> if $ML_y1 == 0
>
>         quietly replace `y' =
> ln[1/sqrt(2*_pi*Sigma_hat^2)]-(1/(2*Sigma_hat^2))*(`y'+Alpha_2-
> Beta_hat*`x')^2
> if $ML_y1 > 0
>
> end
>
>
> I think this will do it. Note that there are some fundamental
> conceptual flaws in what Ben has written above. This can be viewed as
> a linear-form (lf) model of a rather strange sort (given that there
> are two parameter vectors that share an element). I think that could
> be handled by applying constraints that the two slope coefficients
> should be the same.
>
> --------mymlprog.ado--------
> program mymlprog
>         version 10.0
>         args lnf mu1 mu2
>
>         quietly replace `lnf' = log(1/sqrt(2*_pi*$Sigma_hat^2))- ///
>         (1/(2*$Sigma_hat^2))*($ML_y1+`mu1')^2 if $ML_y1 <0
>
>         quietly replace `lnf' =log(normal((`mu2')/$Sigma_hat)- ///
>         normal((`mu1')/$Sigma_hat)) if $ML_y1 == 0
>
>         quietly replace `lnf' =log(1/sqrt(2*_pi*$Sigma_hat^2))- ///
>         (1/(2*$Sigma_hat^2))*($ML_y1+`mu2')^2 ///
>      if $ML_y1 > 0
>
> end
> ---------end of program------
>
> to run, you must first set the global value of Sigma_hat to some
> value (as that does not seem to be a parameter to be estimated)
>
> global Sigma_hat 10
>
> I created a dataset with the desired characteristics with
>
> webuse auto
> g mpgzap = mpg - 25
> ml model lf mymlprog (mpgzap = price) (mpgzap=price)
> ml check
> ml maximize
>
>
>                                                    Number of obs
> =         74
>                                                    Wald chi2(1)
> =       4.76
> Log likelihood = -243.23578                       Prob > chi2
> =     0.0292
>
> ------------------------------------------------------------------------
> ------
>               |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
> Interval]
> -------------
> +----------------------------------------------------------------
> eq1          |
>         price |   .0008808   .0004039     2.18   0.029     .
> 0000891    .0016725
>         _cons |  -2.170081    2.78924    -0.78   0.437
> -7.63689    3.296728
> -------------
> +----------------------------------------------------------------
> eq2          |
>         price |   .0050802   .0023531     2.16   0.031     .
> 0004683    .0096922
>         _cons |  -19.44901   10.88796    -1.79   0.074
> -40.78903    1.891009
> ------------------------------------------------------------------------
> ------
>
> If we now do
>
> cons 1 [eq1]price = [eq2]price
> ml model lf mymlprog (mpgzap = price) (mpgzap=price) , constraints(1)
> ml maximize
>
> the two slope parameters will be equated;
>
>                                                    Number of obs
> =         74
>                                                    Wald chi2(0)
> =          .
> Log likelihood = -246.37235                       Prob > chi2
> =          .
>
>   ( 1)  [eq1]price - [eq2]price = 0
> ------------------------------------------------------------------------
> ------
>               |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
> Interval]
> -------------
> +----------------------------------------------------------------
> eq1          |
>         price |   .0010458   .0004001     2.61   0.009     .
> 0002615    .0018301
>         _cons |  -3.476855   2.776426    -1.25   0.210
> -8.918549    1.964839
> -------------
> +----------------------------------------------------------------
> eq2          |
>         price |   .0010458   .0004001     2.61   0.009     .
> 0002615    .0018301
>         _cons |   -.181225     2.8025    -0.06   0.948
> -5.674024    5.311574
> ------------------------------------------------------------------------
> ------
>
> I think this is estimating the model that Ben has specified. As he
> has an independent
> estimate of the model via Excel Solver, it should be easy to check
> that out.
>
> Kit
>
> Kit Baum, Boston College Economics and DIW Berlin
> http://ideas.repec.org/e/pba1.html
> An Introduction to Modern Econometrics Using Stata:
> http://www.stata-press.com/books/imeus.html
>
>
> *
> *   For searches and help try:
> *   http://www.stata.com/support/faqs/res/findit.html
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2021 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index