different results with nlsur quaids model and QUAIDS -Poi's command (2012)

[David Achiles]

Dear Statalist,

I am trying to built my own codes following old Poi's command "nlsur quaids". (2008)

http://www.stata-journal.com/sjpdf.html?articlenum=st0029_1

At the time I am trying to compare my results with new Poi's command "QUAIDS" (2012)

http://www.stata-journal.com/sjpdf.html?articlenum=st0268

I am trying to compare outcomes step by step with adding each variables one by one. With not quadratic function

replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + `g13'*`lnp3' + `g14'*`lnp4' +  `b1'*(`lnm' - `lnpindex')

the results are similar while I am adding quadratic function,[CODE]replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + `g13'*`lnp3' + `g14'*`lnp4' + `b1'*(`lnm' - `lnpindex') + `l1'/`bofp'*(`lnm' - `lnpindex')^2)[/CODE] the results come with differences coefficients.

My purpose with issue of comparing is that to be sure I am coding on right way.

I highly appreciate your help or recommendation in this issue, might be my way of adding the quadratic terms is wrong or.....??

Thank you


I have STATA 13.1

I am using example data from

webuse food

In my codes, I included also demographics characteristics

My codes:

   clear all
webuse food, clear
   
  *** create dummy variables randomly ****
generate x1 = int(runiform()*4)
generate x2 = (runiform() > 0.7)
 
  rename lnexp lnm
  
**********************************************
***********************MODEL *******************
*****************************************************
set trace off
program nlsurquaids
version 13.1
syntax varlist(min=10 max=10) if, at(name)
tokenize `varlist'
args w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnm x1 x2
   
  tempname a1 a2 a3 a4
scalar `a1' = `at'[1,1]
scalar `a2' = `at'[1,2]
scalar `a3' = `at'[1,3]
scalar `a4' =  1-`a1'-`a2'-`a3'
***************************************************
tempname b1 b2 b3 b4
scalar `b1' = `at'[1,4]
scalar `b2' = `at'[1,5]
scalar `b3' = `at'[1,6]
scalar `b4' = -`b1'-`b2'-`b3'
***************************************************
tempname g11 g12 g13 g14
tempname g21 g22 g23 g24
tempname g31 g32 g33 g34
tempname g41 g42 g43 g44
   
  **************************************************
scalar `g11' = `at'[1,7]
scalar `g12' = `at'[1,8]
scalar `g13' = `at'[1,9]
scalar `g14' = -`g11'-`g12'-`g13'
*************************************************
scalar `g21' = `g12'
scalar `g22' = `at'[1,10]
scalar `g23' = `at'[1,11]
scalar `g24' = -`g21'-`g22'-`g23'
   
  *************************************************
scalar `g31' = `g13'
scalar `g32' = `g23'
scalar `g33' = `at'[1,12]
scalar `g34' = -`g31'-`g32'-`g33'
   
  *************************************************
scalar `g41' = `g14'
scalar `g42' = `g24'
scalar `g43' = `g34'
scalar `g44' = -`g41'-`g42'-`g43'
   
  *************************************************
tempname l1 l2 l3 l4
scalar `l1' = `at'[1,13]
scalar `l2' = `at'[1,14]
scalar `l3' = `at'[1,15]
scalar `l4' = -`l1'-`l2'-`l3'
**********Household demographics****************
tempname r11 r12
tempname r21 r22
tempname r31 r32
 
scalar `r11' = `at'[1,16]
scalar `r12' = `at'[1,17]
scalar `r21' = `at'[1,18]
scalar `r22' = `at'[1,19]
scalar `r31' = `at'[1,20]
scalar `r32' = `at'[1,21]
   
  ****************************************************************************************************
quietly {
tempvar lnpindex
gen double `lnpindex' = (-0.5 + `a1'*`lnp1' + `a2'*`lnp2' + `a3'*`lnp3' + `a4'*`lnp4')
forvalues i = 1/4 {
forvalues j = 1/4 {
replace `lnpindex' = `lnpindex' + 0.5*`g`i'`j''*`lnp`i''*`lnp`j''
}
}
///The b(p) term in the QUAIDS model:
tempvar bofp
gen double `bofp' = 0
forvalues i = 1/4 {
replace `bofp' = `bofp' + `lnp`i''*`b`i''
}
replace `bofp' = exp(`bofp')
****************************************************************************************************************************************
replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + `g13'*`lnp3' + `g14'*`lnp4' +  ///
`b1'*(`lnm' - `lnpindex') + `l1'/`bofp'*(`lnm' - `lnpindex')^2 + `r11'*`x1' +`r12'*`x2')
replace `w2' = (`a2' + `g21'*`lnp1' + `g22'*`lnp2' + `g23'*`lnp3' + `g24'*`lnp4' +  ///
`b2'*(`lnm' - `lnpindex') + `l2'/`bofp'*(`lnm' - `lnpindex')^2 + `r21'*`x1' +`r22'*`x2')
replace `w3' = (`a3' + `g31'*`lnp1' + `g32'*`lnp2' + `g33'*`lnp3' + `g34'*`lnp4' +  ///
`b3'*(`lnm' - `lnpindex') + `l3'/`bofp'*(`lnm' - `lnpindex')^2 + `r31'*`x1' +`r32'*`x2')
   
  }
end
   
  nlsur quaids @ w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnm x1 x2, ifgnls nequations(3) param(a1 a2 a3 ///
g11 g21 g31 g22 g32 g33 b1 b2 b3 l1 l2 l3 r11 r12 r13 r21 r22 r23)

Thank you

[David Achiles]

Sorry for technical mistake

Here Quadratic function

 replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + `g13'*`lnp3' + `g14'*`lnp4' +  ///
`b1'*(`lnm' - `lnpindex') + `l1'/`bofp'*(`lnm' - `lnpindex')^2