How to solve a system of 8 nonlinear equations?

[Hera Li]

Dear All,

I'm kinda new to STATA and I would like to solve a system of 8 nonlinear equations with 8 unknowns. However, I cannot get a RSS that is smaller than 10^-6 after I randomly tried different initial values. Can you give me some suggestions? I really appreciate it!! Below is my code:

program define nlsolvepar

syntax varlist(min=1 max=1) [if], at(name) ///
pi1(real) pi2(real) pi3(real) pi4(real) pi5(real) pi6(real) pi7(real) pi8(real) pu(real)

tempname beta0 beta1 beta2 beta3 alpha0 alpha1 alpha2 alpha3
scalar `beta0' = `at'[1, 1]
scalar `beta1' = `at'[1, 2]
scalar `beta2' = `at'[1, 3]
scalar `beta3' = `at'[1, 4]
scalar `alpha0' = `at'[1, 5]
scalar `alpha1' = `at'[1, 6]
scalar `alpha2' = `at'[1, 7]
scalar `alpha3' = `at'[1, 8]

tempvar yh

* pi1
gen double `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi1'+1 in 1
* pi2
replace `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi2' in 2
*pi3
replace `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`alpha0'+2*`alpha1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi3' in 3
*pi4
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi4' in 4
*pi5
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi5' in 5
*pi6
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`alpha0'+2*`alpha1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi6' in 6
*pi7
replace `yh' = exp(`alpha0'+`alpha1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi7' in 7
*pi8
replace `yh' = exp(`alpha0'+`alpha1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi8' in 8

replace `varlist' = `yh'

end


clear
set obs 8
generate y = 0
replace y = 1 in 1

nl solvepar @ y, pi1(0.1) pi2(0.000001) pi3(0.000001) pi4(0.000001) pi5(0.1) pi6(0.7) pi7(0.000001) pi8(0.000001) pu(0.5) ///
parameters(beta0 beta1 beta2 beta3 alpha0 alpha1 alpha2 alpha3) initial(beta0 3 beta1 3 ///
beta2 3 beta3 3 alpha0 3 alpha1 3 alpha2 3 alpha3 3)

Thanks!
Hera

[William Lisowski]

Welcome to Statalist.

In the output of help nl there is a section describing optimization options. Perhaps one of those will be useful.

[Hera Li]

Thank you for the response! I think the default optimization options work for me. The problem I'm encountering now is that I cannot get a RSS that is smaller enough. The smallest I got was 0.01 although I set different initial values. How can I get a RSS that is small enough?

[William Lisowski]

I think the default optimization options work for me.

I think your results suggest they do not. From help nl

eps(#) specifies the convergence criterion for successive parameter estimates and for the residual sum of squares. The default is eps(1e-5).

[Hera Li]

Yes, the default is 1e-5 and I do want my RSS to be smaller than e-5. However, my problem is that the smallest RSS it can reach is 0.01 and then the iterations stopped. It cannot go any further to reach e-5.

[William Lisowski]

I misunderstood post #1

I cannot get a RSS that is smaller than 10^-6

to mean that 10^-6 was the smallest RSS you could get, and you wanted something smaller (which I thought was perhaps overambitious). I overlooked in post #5 that you said the smallest RSS it can reach is 0.01.

Reviewing your code, I see four occurrences of things like

`p u'
`alp ha0'

which of course is not what you intended. If that is not what you have in your do-file, then there were problems posting your code to Statalist, and we cannot reproduce your results.

Please take a few moments to review the Statalist FAQ linked to from the top of the page, as well as from the Advice on Posting link on the page you used to create your post. Note especially sections 9-12 on how to best pose your question. It's particularly helpful to copy code from your Stata Do-file Editor window, and commands and output from your Stata Results window, and paste them into your Statalist post using code delimiters [CODE] and [/CODE], as described in section 12 of the FAQ.

[Hera Li]

I did copy my codes directly from my do file and there are no such blanks in my codes. The blanks were all presented in the same place, so I assumed it might be from the webpage. Anyways, I removed the spaces.

program define nlsolvepar

syntax varlist(min=1 max=1) [if], at(name) ///
pi1(real) pi2(real) pi3(real) pi4(real) pi5(real) pi6(real) pi7(real) pi8(real) pu(real)

tempname beta0 beta1 beta2 beta3 alpha0 alpha1 alpha2 alpha3
scalar `beta0' = `at'[1, 1]
scalar `beta1' = `at'[1, 2]
scalar `beta2' = `at'[1, 3]
scalar `beta3' = `at'[1, 4]
scalar `alpha0' = `at'[1, 5]
scalar `alpha1' = `at'[1, 6]
scalar `alpha2' = `at'[1, 7]
scalar `alpha3' = `at'[1, 8]

tempvar yh

* pi1
gen double `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi1'+1 in 1
* pi2
replace `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi2' in 2
*pi3
replace `yh' = 1/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`alpha0'+2*`alpha1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
1/(exp(`beta0'+`beta1'+`beta2'+`beta3')+exp(`alpha0' +`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi3' in 3
*pi4
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi4' in 4
*pi5
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi5' in 5
*pi6
replace `yh' = exp(`beta0'+`beta1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`alpha0'+2*`alpha1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`beta0'+`beta1'+`beta2'+`beta3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi6' in 6
*pi7
replace `yh' = exp(`alpha0'+`alpha1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
1/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
1/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+exp(`alp ha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`pu'-`pi7' in 7
*pi8
replace `yh' = exp(`alpha0'+`alpha1')/(exp(`beta0'+`beta1')+exp(`alpha0'+`alpha1')+1)* ///
exp(`beta0'+2*`beta1')/(exp(`beta0'+2*`beta1')+exp(`alpha0'+2*`alpha1')+1 )*(1-`pu')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')/(exp(`beta0'+`beta1'+`beta2'+`beta3')+ ///
exp(`alpha0'+`alpha1'+`alpha2'+`alpha3')+1)* ///
exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')/(exp(`beta0'+2*`beta1'+`beta2'+2*`beta3')+ ///
exp(`alpha0'+2*`alpha1'+`alpha2'+2*`alpha3')+1)*`p u'-`pi8' in 8

replace `varlist' = `yh'

end


clear
set obs 8
generate y = 0
replace y = 1 in 1

nl solvepar @ y, pi1(0.1) pi2(0.000001) pi3(0.000001) pi4(0.000001) pi5(0.1) pi6(0.7) pi7(0.000001) pi8(0.000001) pu(0.5) ///
parameters(beta0 beta1 beta2 beta3 alpha0 alpha1 alpha2 alpha3) initial(beta0 3 beta1 3 ///
beta2 3 beta3 3 alpha0 3 alpha1 3 alpha2 3 alpha3 3)

[Hera Li]

And when I repost it, the same problem occurred again. I guess something wrong with the uploading.

[William Lisowski]

To assure maximum readability of code that you post, which I hope will also avoid the problems with your posting, please copy your code from the Do-file editor, or whatever editor you use, and paste it into a code block in the Forum editor using code delimiters [CODE] and [/CODE], as explained in section 12 of the Statalist FAQ linked to at the top of the page that I referenced in post #6. For example, the following:

[CODE]
. sysuse auto, clear
(1978 Automobile Data)

. describe make price

storage display value
variable name type format label variable label
-----------------------------------------------------------------
make str18 %-18s Make and Model
price int %8.0gc Price
[/CODE]

will be presented in the post as the following:

. sysuse auto, clear
(1978 Automobile Data)

. describe make price

              storage   display    value
variable name   type    format     label      variable label
-----------------------------------------------------------------
make            str18   %-18s                 Make and Model
price           int     %8.0gc                Price