Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: The another problem about mm_root


From   Matthew Baker <matthew.baker@hunter.cuny.edu>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: The another problem about mm_root
Date   Tue, 6 Sep 2011 09:37:19 -0400

Without knowing what the data in the v vector look like, it is
impossible to tell exactly what is happening, but it looks like you
might be having one or two mathematical problems, not problems with
mm_root. The given function is undefined at x=0. This is due to the
very last component of the function,  b11*(f12-x)/(x*(1+x)^11), where
x appears in the denominator.

The help for mm_root says that if the function is undefined at some
point in the range, mm_root will return x=. If I run the following
code, with the range modified to be strictly positive (and with
strictly random data):

/* Begin example */
clear all
mata
set seed 8675309
v=invnormal(runiform(100,26))
function y(x,v) {
p=v[2];b0=v[3];b1=v[4];b2=v[5];b3=v[6];b4=v[7];b5=v[8];b6=v[9];b7=v[10];b8=v[11];b9=v[12];b10=v[13];b11=v[14];f1=v[15];f2=v[16];f3=v[17];f4=v[18];f5=v[19];f6=v[20];f7=v[21];f8=v[22];f9=v[23];f10=v[24];f11=v[25];f12=v[26]
return(-p+b0+b0*(f1-x)/(1+x)+b1*(f2-x)/(1+x)^2+b2*(f3-x)/(1+x)^3+b3*(f4-x)/(1+x)^4+b4*(f5-x)/(1+x)^5+b5*(f6-x)/(1+x)^6+b6*(f7-x)/(1+x)^7+b7*(f8-x)/(1+x)^8+b8*(f9-x)/(1+x)^9+b9*(f10-x)/(1+x)^10+b10*(f11-x)/(1+x)^11+b11*(f12-x)/(x*(1+x)^11))
}
Z=J(rows(v),2,.)
for (i=1;i<=rows(v);i++) {
r=mm_root(x=.,&y(),0.01,1,1e-9,1000,v[i,.])
Z[i,1]=x
Z[i,2]=r
}
end
/* end example */

I get roots to the given function about half of the time (mm_root
returns 0), while the other half of the time the function doesn't
cross the zero line (mm_root returns 3). So, I think the key problem
is the inclusion of zero in the range, but it could also be that the
data does not allow the function to cross the origin.

Hope that helps!

Matt Baker




-- 
Dr. Matthew J. Baker
Department of Economics
Hunter College and the Graduate Center, CUNY

As a simple experiment, consider

On Sat, Sep 3, 2011 at 9:16 AM, dai_yunhao <dai_yunhao@smail.hust.edu.cn> wrote:
> Dear Everyone,
>
> Thank you for your suggestion, and I try to use vector to instead for the
> argument, there is no error indeed, but this time I face another problem.
> My program as follows:
>
> clear all
> use gls.dta,clear
> gen x=.
> mata
> v=J(1,1,.)
> st_view(v,., "x price b0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 froe_1 froe_2
> froe_3 froe_4 froe_5 froe_6 froe_7 froe_8 froe_9 froe_10 froe_11 froe_12")
> function y(x,v) {
> p=v[2];b0=v[3];b1=v[4];b2=v[5];b3=v[6];b4=v[7];b5=v[8];b6=v[9];b7=v[10];b8=v[11];b9=v[12];b10=v[13];b11=v[14];f1=v[15];f2=v[16];f3=v[17];f4=v[18];f5=v[19];f6=v[20];f7=v[21];f8=v[22];f9=v[23];f10=v[24];f11=v[25];f12=v[26]
> return(-p+b0+b0*(f1-x)/(1+x)+b1*(f2-x)/(1+x)^2+b2*(f3-x)/(1+x)^3+b3*(f4-x)/(1+x)^4+b4*(f5-x)/(1+x)^5+b5*(f6-x)/(1+x)^6+b6*(f7-x)/(1+x)^7+b7*(f8-x)/(1+x)^8+b8*(f9-x)/(1+x)^9+b9*(f10-x)/(1+x)^10+b10*(f11-x)/(1+x)^11+b11*(f12-x)/(x*(1+x)^11))
> }
> for (i=1;i<=rows(v);i++) {
> r=mm_root(x=.,&y(),0,1,1e-9,1000,v[i,.])
> v[i,1]=x
> }
> end
>
> When I run the program, the value of x is always the missing value ".", I
> try to change the interval
> e.g.
> r=mm_root(x=.,&y(),*-10,10*,1e-9,1000,v[i,.])
>
> but the value of x is -10 or 10, apparently it is not the real valut of x.
> Could anyone tell me what's wrong with the program?
> Thanks for your help!
>
> Dai Yunhao
>
>
> --
> View this message in context: http://statalist.1588530.n2.nabble.com/The-another-problem-about-mm-root-tp6756596p6756596.html
> Sent from the Statalist mailing list archive at Nabble.com.
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index