*********************************************************************

* Replication file 10/11 for "Are Sports Teams Multiproduct Firms?"
* by K.G. Stewart and J.C.H. Jones (2010) Empirical Economics 39(2),
* 487-514.

* Estimation of MSGM factor demand systems for baseball player
* characteristics, the model M0. We estimate the nonseparable model
* with concavity imposed.
* This printout compute marginal costs and use them to obtain the 
* product-specific returns to scale values reported in Table 8, 
* including their standard errors.

********************************************************************* 

* Sample consists of 26 teams x 6 seasons (1986-91) 
SAMPLE 1 156;

* Read in team characteristics; documented in Team.doc
read(Team.dat) season team league attend attlag ticprice &
pmarkup gameswon gameslst wins salaries nplayers strkouts walks &
slugavg nhitters npitchrs experhit experpit nstars roadatt &
aleast alcentrl alwest nleast nlcentrl nlwest
* Delete variables not used in this analysis
delete season team league attlag ticprice pmarkup gameswon gameslst salaries nplayers

* read in hedonic prices
read(HedPrice1.dat) p1 p2 p3 p4 nobs

* define team outputs
* This scaling of road attendance yields fastest convergence
GENR Y1=WINS
GENR Y2=ROADatt/10000000

* conference dummies
genr dum1=aleast
genr dum2=alcentrl
genr dum3=alwest
genr dum4=nleast
genr dum5=nlcentrl  

* define team factor inputs: experience, hitting, pitching, stars
GENR q1=experhit+experpit
GENR q2=nhitters*slugavg
GENR q3=npitchrs*strkouts/walks
GENR q4=nstars

* DEFLATE HEDONIC PRICES TO 1991 DOLLARS
SAMPLE 1 26      
GENR p1=p1*135.0/109.6
GENR p2=p2*135.0/109.6
GENR p3=p3*135.0/109.6
GENR p4=p4*135.0/109.6
SAMPLE 27 52
GENR p1=p1*135.0/113.6
GENR p2=p2*135.0/113.6
GENR p3=p3*135.0/113.6
GENR p4=p4*135.0/113.6
SAMPLE 53 78
GENR p1=p1*135.0/118.3
GENR p2=p2*135.0/118.3
GENR p3=p3*135.0/118.3
GENR p4=p4*135.0/118.3
SAMPLE 79 104
GENR p1=p1*135.0/124.0
GENR p2=p2*135.0/124.0
GENR p3=p3*135.0/124.0
GENR p4=p4*135.0/124.0
SAMPLE 105 130
GENR p1=p1*135.0/130.7
GENR p2=p2*135.0/130.7
GENR p3=p3*135.0/130.7
GENR p4=p4*135.0/130.7

SAMPLE 1 156

* Generate cost identity and cost shares 
* (This synthetic (hedonic) cost series is used 
* only for descriptive purposes in calculating the implied cost shares, not
* in estimation of cost function parameters. It can be compared with fitted
* cost series predicted by estimated model.)
genr cost=q1*p1+q2*p2+q3*p3+q4*p4 
genr s1=q1*p1/cost                
genr s2=q2*p2/cost                
genr s3=q3*p3/cost                
genr s4=q4*p4/cost

stat Y1 Y2 p1 p2 p3 p4 s1 s2 s3 s4

* Obtain mean factor inputs:
* qm1=135.73077
* qm2=5.35117
* qm3=17.92278
* qm4=2.61538
stat q1 / mean=qm1
stat q2 / mean=qm2
stat q3 / mean=qm3
stat q4 / mean=qm4

GENR DENOM=qm1*p1+qm2*p2+qm3*p3+qm4*p4

gen1 theta1=qm1/2
gen1 theta2=qm2/2
gen1 theta3=qm3/2
gen1 theta4=qm4/2
               
* generate price constructs used in estimation

GENR Pd1=p1/DENOM
GENR Pd2=p2/DENOM
GENR Pd3=p3/DENOM
GENR Pd4=p4/DENOM

GENR Pd11=pd1*pd1
GENR Pd12=pd1*pd2
GENR Pd13=pd1*pd3
GENR Pd14=pd1*pd4
GENR Pd21=pd2*pd1
GENR Pd22=pd2*pd2
GENR Pd23=pd2*pd3
GENR Pd24=pd2*pd4
GENR Pd31=pd3*pd1
GENR Pd32=pd3*pd2
GENR Pd33=pd3*pd3
GENR Pd34=pd3*pd4
GENR Pd41=pd4*pd1
GENR Pd42=pd4*pd2
GENR Pd43=pd4*pd3
GENR Pd44=pd4*pd4

GENR Pd11=Pd11-Pd14-(Pd41-Pd44)
GENR Pd12=Pd12-Pd14-(Pd42-Pd44)
GENR Pd13=Pd13-Pd14-(Pd43-Pd44)
GENR Pd21=Pd21-Pd24-(Pd41-Pd44)
GENR Pd22=Pd22-Pd24-(Pd42-Pd44)
GENR Pd23=Pd23-Pd24-(Pd43-Pd44)
GENR Pd31=Pd31-Pd34-(Pd41-Pd44)
GENR Pd32=Pd32-Pd34-(Pd42-Pd44)
GENR Pd33=Pd33-Pd34-(Pd43-Pd44)

GENR Pd12=2*Pd12
GENR Pd13=2*Pd13
GENR Pd23=2*Pd23

GENR Pd1=Pd1-Pd4
GENR Pd2=Pd2-Pd4
GENR Pd3=Pd3-Pd4

* Obtain parallel scalar means used in reconstruction of 
* total cost predicted at pt of variable means

stat p1 / mean=m1
stat p2 / mean=m2
stat p3 / mean=m3
stat p4 / mean=m4

GEN1 m11=m1*m1
GEN1 m12=m1*m2 
GEN1 m13=m1*m3 
GEN1 m14=m1*m4 
GEN1 m21=m2*m1 
GEN1 m22=m2*m2 
GEN1 m23=m2*m3 
GEN1 m24=m2*m4 
GEN1 m31=m3*m1 
GEN1 m32=m3*m2 
GEN1 m33=m3*m3 
GEN1 m34=m3*m4 
GEN1 m41=m4*m1 
GEN1 m42=m4*m2 
GEN1 m43=m4*m3 
GEN1 m44=m4*m4 
                         
GEN1 m11=m11-m14-(m41-m44)
GEN1 m12=m12-m14-(m42-m44)
GEN1 m13=m13-m14-(m43-m44)
GEN1 m21=m21-m24-(m41-m44)
GEN1 m22=m22-m24-(m42-m44)
GEN1 m23=m23-m24-(m43-m44)
GEN1 m31=m31-m34-(m41-m44)
GEN1 m32=m32-m34-(m42-m44)
GEN1 m33=m33-m34-(m43-m44) 

* check symmetry: pairs should be equal
print m12 m21 
print m13 m31
print m23 m32 

GEN1 m12=2*m12
GEN1 m13=2*m13
GEN1 m23=2*m23

GEN1 mDENOM=qm1*m1+qm2*m2+qm3*m3+qm4*m4
GEN1 mdenom2=mdenom*mdenom

GEN1 md1=m1/mDENOM
GEN1 md2=m2/mDENOM
GEN1 md3=m3/mDENOM
GEN1 md4=m4/mDENOM

* Obtain parallel scalar means used in reconstruction of 
* factor demands and shares predicted at pt of variable means

GEN1 md11=md1*md1
GEN1 md12=md1*md2 
GEN1 md13=md1*md3 
GEN1 md14=md1*md4 
GEN1 md21=md2*md1 
GEN1 md22=md2*md2 
GEN1 md23=md2*md3 
GEN1 md24=md2*md4 
GEN1 md31=md3*md1 
GEN1 md32=md3*md2 
GEN1 md33=md3*md3 
GEN1 md34=md3*md4 
GEN1 md41=md4*md1 
GEN1 md42=md4*md2 
GEN1 md43=md4*md3 
GEN1 md44=md4*md4 
                         
GEN1 md11=md11-md14-(md41-md44)
GEN1 md12=md12-md14-(md42-md44)
GEN1 md13=md13-md14-(md43-md44)
GEN1 md21=md21-md24-(md41-md44)
GEN1 md22=md22-md24-(md42-md44)
GEN1 md23=md23-md24-(md43-md44)
GEN1 md31=md31-md34-(md41-md44)
GEN1 md32=md32-md34-(md42-md44)
GEN1 md33=md33-md34-(md43-md44)  

GEN1 md12=2*md12
GEN1 md13=2*md13
GEN1 md23=2*md23

gen1 md1=md1-md4
gen1 md2=md2-md4
gen1 md3=md3-md4

stat y1 / mean=ym1
stat y2 / mean=ym2

** INITIAL PARAMETERIZATION **
* Shazam requires that the theta's (/2) appear numerically
* Multiplicative factor 2 is included in parameter D1

nl 4 / ncoef=41 iter=999 
Eq q1=(S11*Pd1+S12*Pd2+S13*Pd3            &
     -67.8654*(S11*Pd11+S12*Pd12+S13*Pd13 &
                       +S22*Pd22+S23*Pd23 &
                                +S33*Pd33))*(Y1+BETA*Y2) &
     +d11*dum1+d12*dum2+d13*dum3+d14*dum4+d15*dum5 &
     +B11*(Y1+BETA*Y2)+B1+LAMBDA1*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q2=(S12*Pd1+S22*Pd2+S23*Pd3            &
     -2.67558*(S11*Pd11+S12*Pd12+S13*Pd13 &
                       +S22*Pd22+S23*Pd23 &
                                +S33*Pd33))*(Y1+BETA*Y2) &
     +d21*dum1+d22*dum2+d23*dum3+d24*dum4+d25*dum5 &
     +B22*(Y1+BETA*Y2)+B2+LAMBDA2*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q3=(S13*Pd1+S23*Pd2+S33*Pd3            &
     -8.96111*(S11*Pd11+S12*Pd12+S13*Pd13 &
                       +S22*Pd22+S23*Pd23 &
                                +S33*Pd33))*(Y1+BETA*Y2) &
     +d31*dum1+d32*dum2+d33*dum3+d34*dum4+d35*dum5 &
     +B33*(Y1+BETA*Y2)+B3+LAMBDA3*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q4=-((S11+S12+S13)*Pd1+(S12+S22+S23)*Pd2+(S13+S23+S33)*Pd3  &
     +1.30769*(S11*Pd11+S12*Pd12+S13*Pd13 &
                       +S22*Pd22+S23*Pd23 &
                                +S33*Pd33))*(Y1+BETA*Y2) &
     +d41*dum1+d42*dum2+d43*dum3+d44*dum4+d45*dum5 &
     +B44*(Y1+BETA*Y2)+B4+LAMBDA4*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
end

* Check integrity of components of the cost expression that will be used
* in elasticities by computing:
* Cost function at pt of variable means
test ((S11*m11+S12*m12+S13*m13 &
              +S22*m22+S23*m23 &
                      +S33*m33)/(2*mdenom) &
     +b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2)+b1*m1+b2*m2+b3*m3+b4*m4 &
     +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2)
   
* pscale1=1-((C(y1,y2)-C(0,y2))/(MC1*y1))
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2) &
      +(b1*m1+b2*m2+b3*m3+b4*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2) &
      -((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(beta*ym2) &
      -(b1*m1+b2*m2+b3*m3+b4*m4) &
      -(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D2*Ym2*Ym2) )*ym1/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(2*Ym1+D1*Ym2) )
* pscale1: sensitivity to reparameterization
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*ym1 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2) )*ym1/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(2*Ym1+D1*Ym2) )

* pscale2=1-((C(y1,y2)-C(y1,0))/(MC2*y2))
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2) &
      +(b1*m1+b2*m2+b3*m3+b4*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2) &
      -((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*ym1 &
      -(b1*m1+b2*m2+b3*m3+b4*m4) &
      -(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1) )*ym2/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1+2*D2*Ym2) )
* pscale2: sensitivity to reparameterization
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta*ym2 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1*Ym2+D2*Ym2*Ym2) )*ym2/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1+2*D2*Ym2) )

** ALTERNATIVE (BUT EQUIVALENT) PARAMETERIZATION **        

GENR Z11=Pd1-THETA1*Pd11
GENR Z12=Pd2-THETA1*Pd12
GENR Z13=Pd3-THETA1*Pd13
GENR Z14=-THETA1*Pd22
GENR Z15=-THETA1*Pd23
GENR Z16=-THETA1*Pd33
GENR Z21=-THETA2*Pd11
GENR Z22=Pd1-THETA2*Pd12
GENR Z23=-THETA2*Pd13
GENR Z24=Pd2-THETA2*Pd22
GENR Z25=Pd3-THETA2*Pd23
GENR Z26=-THETA2*Pd33
GENR Z31=-THETA3*Pd11
GENR Z32=-THETA3*Pd12
GENR Z33=Pd1-THETA3*Pd13
GENR Z34=-THETA3*Pd22
GENR Z35=Pd2-THETA3*Pd23
GENR Z36=Pd3-THETA3*Pd33
GENR Z41=-(Pd1+THETA4*Pd11)
GENR Z42=-(Pd1+Pd2+THETA4*Pd12)
GENR Z43=-(Pd1+Pd3+THETA4*Pd13)
GENR Z44=-(Pd2+THETA4*Pd22)
GENR Z45=-(Pd2+Pd3+THETA4*Pd23)
GENR Z46=-(Pd3+THETA4*Pd33)

* Multiplicative factor 2 is included in parameter D1  
nl 4 / ncoef=41 iter=999
Eq q1=(Y1+BETA*Y2)                                               &
          *(B11+S11*Z11+S12*Z12+S13*Z13+S22*Z14+S23*Z15+S33*Z16) &
          +d11*dum1+d12*dum2+d13*dum3+d14*dum4+d15*dum5 &
          +B1+LAMBDA1*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q2=(Y1+BETA*Y2)                                               &
          *(B22+S11*Z21+S12*Z22+S13*Z23+S22*Z24+S23*Z25+S33*Z26) &
          +d21*dum1+d22*dum2+d23*dum3+d24*dum4+d25*dum5 &
          +B2+LAMBDA2*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q3=(Y1+BETA*Y2)                                               &
          *(B33+S11*Z31+S12*Z32+S13*Z33+S22*Z34+S23*Z35+S33*Z36) &
          +d31*dum1+d32*dum2+d33*dum3+d34*dum4+d35*dum5 &
          +B3+LAMBDA3*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q4=(Y1+BETA*Y2)                                               &
          *(B44+S11*Z41+S12*Z42+S13*Z43+S22*Z44+S23*Z45+S33*Z46) &
          +d41*dum1+d42*dum2+d43*dum3+d44*dum4+d45*dum5 &
          +B4+LAMBDA4*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
end

* Check integrity of components of the cost expression that will be used
* in elasticities by computing:
* Cost function at pt of variable means
test ((S11*m11+S12*m12+S13*m13 &
              +S22*m22+S23*m23 &
                      +S33*m33)/(2*mdenom) &
     +b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2)+b1*m1+b2*m2+b3*m3+b4*m4 &
     +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2)
     
* pscale1=1-((C(y1,y2)-C(0,y2))/(MC1*y1))
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2) &
      +(b1*m1+b2*m2+b3*m3+b4*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2) &
      -((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(beta*ym2) &
      -(b1*m1+b2*m2+b3*m3+b4*m4) &
      -(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D2*Ym2*Ym2) )*ym1/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(2*Ym1+D1*Ym2) )
* pscale1: sensitivity to reparameterization
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*ym1 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2) )*ym1/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(2*Ym1+D1*Ym2) )

* pscale2=1-((C(y1,y2)-C(y1,0))/(MC2*y2))
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*(ym1+beta*ym2) &
      +(b1*m1+b2*m2+b3*m3+b4*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2+D2*Ym2*Ym2) &
      -((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*ym1 &
      -(b1*m1+b2*m2+b3*m3+b4*m4) &
      -(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1) )*ym2/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1+2*D2*Ym2) )
* pscale2: sensitivity to reparameterization
test 1- ( ((S11*m11+S12*m12+S13*m13 &
                   +S22*m22+S23*m23 &
                           +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta*ym2 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1*Ym2+D2*Ym2*Ym2) )*ym2/ &
     ( ((S11*m11+S12*m12+S13*m13 &
                +S22*m22+S23*m23 &
                        +S33*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1+2*D2*Ym2) )

** CONCAVITY IMPOSED USING WILEY-SCHMIDT-BRAMBLE REPARAMETERIZATION **
* Multiplicative factor 2 is included in parameter D1: this turns out to be
* critical to the quick convergence of the concavity-constrained system   
nl 4 / ncoef=41 iter=999 coef=g predict=qpred
Eq q1=-(G11*G11*Z11+G11*G21*Z12+G11*G31*Z13+(G21*G21+G22*G22)*Z14 &
         +(G31*G21+G32*G22)*Z15+(G31*G31+G32*G32+G33*G33)*Z16-B11)*(Y1+BETA*Y2) &
      +d11*dum1+d12*dum2+d13*dum3+d14*dum4+d15*dum5 &   
      +B1+LAMBDA1*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q2=-(G11*G11*Z21+G11*G21*Z22+G11*G31*Z23+(G21*G21+G22*G22)*Z24 &
         +(G31*G21+G32*G22)*Z25+(G31*G31+G32*G32+G33*G33)*Z26-B22)*(Y1+BETA*Y2) &
      +d21*dum1+d22*dum2+d23*dum3+d24*dum4+d25*dum5 &   
      +B2+LAMBDA2*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q3=-(G11*G11*Z31+G11*G21*Z32+G11*G31*Z33+(G21*G21+G22*G22)*Z34 &
         +(G31*G21+G32*G22)*Z35+(G31*G31+G32*G32+G33*G33)*Z36-B33)*(Y1+BETA*Y2) &
      +d31*dum1+d32*dum2+d33*dum3+d34*dum4+d35*dum5 &   
      +B3+LAMBDA3*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
Eq q4=-(G11*G11*Z41+G11*G21*Z42+G11*G31*Z43+(G21*G21+G22*G22)*Z44 &
         +(G31*G21+G32*G22)*Z45+(G31*G31+G32*G32+G33*G33)*Z46-B44)*(Y1+BETA*Y2) &
      +d41*dum1+d42*dum2+d43*dum3+d44*dum4+d45*dum5 &   
      +B4+LAMBDA4*(Y1*Y1+D1*Y1*Y2+D2*Y2*Y2)
end

* We use the following parameter equivalencies
*  s11=-g11*g11                  
*  s12=-g11*g21                  
*  s13=-g11*g31                  
*  s22=-(g21*g21+g22*g22)        
*  s23=-(g21*g31+g22*g32)        
*  s33=-(g31*g31+g32*g32+g33*g33)

* pscale1 (reparameterized)
test 1- ( ((-g11*g11*m11-g11*g21*m12-g11*g31*m13 &          
    -(g21*g21+g22*g22)*m22-(g21*g31+g22*g32)*m23 &
                  -(g31*g31+g32*g32+g33*g33)*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*ym1 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(Ym1*Ym1+D1*Ym1*Ym2) )*ym1/ &
     ( ((-g11*g11*m11-g11*g21*m12-g11*g31*m13 &          
 -(g21*g21+g22*g22)*m22-(g21*g31+g22*g32)*m23 &
               -(g31*g31+g32*g32+g33*g33)*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4) &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(2*Ym1+D1*Ym2) )
* pscale2 (reparameterized)
test 1- ( ((-g11*g11*m11-g11*g21*m12-g11*g31*m13 &          
    -(g21*g21+g22*g22)*m22-(g21*g31+g22*g32)*m23 &
                  -(g31*g31+g32*g32+g33*g33)*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta*ym2 &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1*Ym2+D2*Ym2*Ym2) )*ym2/ &
        ( ((-g11*g11*m11-g11*g21*m12-g11*g31*m13 &          
    -(g21*g21+g22*g22)*m22-(g21*g31+g22*g32)*m23 &
                  -(g31*g31+g32*g32+g33*g33)*m33)/(2*mdenom)+b11*m1+b22*m2+b33*m3+b44*m4)*beta &
      +(lambda1*m1+lambda2*m2+lambda3*m3+lambda4*m4)*(D1*Ym1+2*D2*Ym2) )      

stop