/****************************************
* EWreg.ado
*
* Runs Error-In-Variable model regression, with 1 mismeasured regressors and several
* perfectly measured regressors.
*
* Uses the High-Order-Moments method of Erickson Whited 2000, Journal of Political
* Economy, also described in Erickson Whited 2002, Econometric Theory.
*
* Author:  Robert Parham, University of Rochester
*          June 30, 2012
*
* Version: 1.0
*
* Based on Gauss code provided by Toni M. Whited
*
*
* USAGE:
* syntax: EWreg depvar misindep [indepvars] [if] [in] [, option]
* options:
*	METHod 		- one of "GMM3","GMM4","GMM5","GMM6","GMM7". default is "GMM3"
*	BXint		- a numlist of starting values for beta (the coefficient on misindep)
*				  zero means use the best guess, and should be part of the list.
* 				  defualt is 0.
*	HAScons 	- indicates that indepvars contains a constant variable, and so a constant
* 				  should not be added in the estimation process
*	NOPRN 		- supress printing of results
*
* RETURN:
* EWreg saves the following in e()
*
* scalars:
* 	e(N) 		- number of observations
*	e(Jstat) 	- J statistic on overidentifying restrictions
*	e(Jval) 	- p-val on Jstat
*	e(dfree)	- Degrees of freedom for Jval
*	e(rho) 		- the rho^2 value (see EW2000 JPE)
*	e(tau) 		- the tau^2 value (see EW2000 JPE)
*	e(SErho) 	- s.e. on rho^2
*	e(SEtau) 	- s.e. on tau^2
* 	e(obj)		- value of GMM objective function
*	e(ID3stat) 	- Test statistic on identification test with 3 moments
*	e(ID3val) 	- p-val on ID3val
*	e(ID4stat) 	- Test statistic on identification test with 4 moments
*	e(ID4val) 	- p-val on ID4val
*
* macros:
*	e(bxint) 	- the numlist of initial guesses for beta
*
* matrices:
* 	e(b) 		- regression coeffiecints
* 	e(V) 		- Variance-Covariance matrix on e(b)
* 	e(inflnc) 	- Influence functions on (misindep,indepvars)
*	e(serr) 	- standard errors on e(b)
* 	e(vcrhotau) - Variance-Covariance matrix on (rho,tau)
*	e(w) 		- weighting matrix used in GMM step of estimation
* 	e(inflncrhotau)	- Influence functions on (rho,tau)
*
****************************************/


capture program drop EWreg
program define EWreg, eclass
	version 10.1
	syntax varlist(min=2 numeric) [if] [in] [, METHod(string) BXint(numlist) HAScons NOPRN]
	marksample touse

	quietly count if `touse'
	if `r(N)' == 0 error 2000

	gettoken depvar varlist: varlist 		// separate varlist
	gettoken misindep indepvars: varlist 	// separate varlist

	if "`bxint'"=="" local bxint = "0"

	tempname b V cst
	mata: doEW("`depvar'", "`misindep'", "`indepvars'", "`method'", "`bxint'", "`hascons'" , "`touse'" )
	mat `b' = r(beta)
	mat `V' = r(VCmat)
	local cst = cond("`hascons'" == "", "_cons", "")
	local vnames `misindep' `cst' `indepvars'
	matname `V' `vnames'
	matname `b' `vnames', c(.)
	local N = r(N)
	ereturn post `b' `V', depname(`depvar') obs(`N') esample(`touse')
	ereturn matrix inflnc  = inflnc
	ereturn matrix serr  = serr
	ereturn matrix vcrhotau = vcrhotau
	ereturn matrix inflncrhotau = inflncrhotau
	ereturn matrix w = w
	ereturn scalar Jstat = r(Jstat)
	ereturn scalar Jval = r(Jval)
	ereturn scalar dfree = r(dfree)
	ereturn scalar rho = r(rho)
	ereturn scalar tau = r(tau)
	ereturn scalar SErho = r(SErho)
	ereturn scalar SEtau = r(SEtau)
	ereturn scalar obj = r(obj)
	ereturn local  bxint = "`bxint'"
	ereturn scalar ID3val = r(ID3val)
	ereturn scalar ID3stat = r(ID3stat)
	ereturn scalar ID4val = r(ID4val)
	ereturn scalar ID4stat = r(ID4stat)
	if ("`noprn'"=="") {
		display _newline "`method' Errors-In-Variables results" _col(60) "Number of obs = " e(N)
		ereturn display
		if ( e(Jstat) > 0 ) {
			display "Sargan-Hansen J statistic: " %7.3f e(Jstat)
			display "Chi-sq(" %3.0f e(dfree) " ) P-val = " ///
			%5.4f e(Jval) _newline
		}
	}
end

version 10.1
mata:
mata clear

// subroutine to define the f vector.
function deff(a,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,ey3_d,ey4x_d,/*
			 */ ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,/*
			 */ ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,eestim,nneeqq)
{
	f=J(1,nneeqq,0)

	f[,1]=ey2_d-(a[1,1]^2)*a[2,1]-a[3,1]
	f[,2]=eyx_d-a[1,1]*a[2,1]
	f[,3]=ex2_d-a[2,1]-a[4,1]
	f[,4]=ey2x_d-(a[1,1]^2)*a[5,1]
	f[,5]=eyx2_d-a[1,1]*a[5,1]

	if (eestim>1)
	{
		f[,6]=ey3x_d-(a[1,1]^3)*a[6,1]-3*a[1,1]*a[2,1]*a[3,1]
		f[,7]=ey2x2_d-(a[1,1]^2)*(a[6,1]+a[2,1]*a[4,1])-a[3,1]*(a[2,1]+a[4,1])
		f[,8]=eyx3_d-a[1,1]*a[6,1]-3*a[1,1]*a[2,1]*a[4,1]
	}

	if (eestim>2)
	{
		f[,9]=ex3_d-a[5,1]-a[7,1]
		f[,10]=ey3_d-(a[1,1]^3)*a[5,1]-a[8,1]
		f[,11]=ey4x_d-(a[1,1]^4)*a[9,1]-6*(a[1,1]^2)*a[5,1]*a[3,1]-4*a[1,1]*a[2,1]*a[8,1]
		f[,12]=ey3x2_d-(a[1,1]^3)*(a[9,1]+a[5,1]*a[4,1])-3*a[1,1]*a[5,1]*a[3,1]-a[8,1]*(a[2,1]+a[4,1])
		f[,13]=ey2x3_d-(a[1,1]^2)*(a[9,1]+3*a[5,1]*a[4,1]+a[2,1]*a[7,1])-a[3,1]*(a[5,1]+a[7,1])
		f[,14]=eyx4_d-a[1,1]*(a[9,1]+6*a[5,1]*a[4,1]+4*a[2,1]*a[7,1])
	}

	if (eestim>3)
	{
		f[,15]=ey4_d-(a[1,1]^4)*a[6,1]-6*(a[1,1]^2)*a[2,1]*a[3,1]-a[10,1]
		f[,16]=ex4_d-a[6,1]-6*a[2,1]*a[4,1]-a[11,1]
		f[,17]=ey5x_d-(a[1,1]^5)*a[12,1]-10*(a[1,1]^3)*a[6,1]*a[3,1]-10*(a[1,1]^2)*a[5,1]*a[8,1]-5*a[1,1]*a[2,1]*a[10,1]
		f[,18]=ey4x2_d-(a[1,1]^4)*(a[12,1]+a[6,1]*a[4,1])-6*(a[1,1]^2)*(a[6,1]*a[3,1]+a[2,1]*a[3,1]*a[4,1])-4*a[1,1]*a[5,1]*a[8,1]-a[10,1]*(a[2,1]+a[4,1])
		f[,19]=ey3x3_d-(a[1,1]^3)*(a[12,1]+3*a[6,1]*a[4,1]+a[5,1]*a[7,1])-a[1,1]*(3*a[6,1]*a[3,1]+9*a[2,1]*a[3,1]*a[4,1])-a[8,1]*(a[5,1]+a[7,1])
		f[,20]=ey2x4_d-(a[1,1]^2)*(a[12,1]+6*a[6,1]*a[4,1]+4*a[5,1]*a[7,1]+a[2,1]*a[11,1])-a[3,1]*(a[6,1]+6*a[2,1]*a[4,1]+a[11,1])
		f[,21]=eyx5_d-a[1,1]*(a[12,1]+10*a[6,1]*a[4,1]+10*a[5,1]*a[7,1]+5*a[2,1]*a[11,1])
	}

	if (eestim>4)
	{
		f[,22]=ey5_d-(a[1,1]^5)*a[9,1]-10*(a[1,1]^3)*a[5,1]*a[3,1]/*
			 */-10*(a[1,1]^2)*a[2,1]*a[8,1]-a[13,1]
		f[,23]=ex5_d-a[9,1]-10*a[5,1]*a[4,1]-10*a[2,1]*a[7,1]-a[14,1]
		f[,24]=ey6x_d-(a[1,1]^6)*a[15,1]/*
			 */-15*(a[1,1]^4)*a[9,1]*a[3,1]/*
			 */-20*(a[1,1]^3)*a[6,1]*a[8,1]/*
			 */-15*(a[1,1]^2)*a[5,1]*a[10,1]/*
			 */-6*a[1,1]*a[2,1]*a[13,1]
		f[,25]=ey5x2_d-(a[1,1]^5)*(a[15,1]+a[9,1]*a[4,1])/*
			 */-10*(a[1,1]^3)*(a[9,1]*a[3,1]+a[5,1]*a[3,1]*a[4,1])/*
			 */-10*(a[1,1]^2)*(a[6,1]*a[8,1]+a[2,1]*a[8,1]*a[4,1])/*
			 */-5*a[1,1]*a[5,1]*a[10,1]/*
			 */-a[13,1]*(a[2,1]+a[4,1])
		f[,26]=ey4x3_d/*
			 */-(a[1,1]^4)*(a[15,1]+3*a[9,1]*a[4,1]+a[6,1]*a[7,1])/*
			 */-6*(a[1,1]^2)*(a[9,1]*a[3,1]+3*a[5,1]*a[3,1]*a[4,1]+a[2,1]*a[3,1]*a[7,1])/*
			 */-4*a[1,1]*(a[6,1]*a[8,1]+3*a[2,1]*a[8,1]*a[4,1])/*
			 */-a[10,1]*(a[5,1]+a[7,1])
		f[,27]=ey3x4_d/*
			 */-(a[1,1]^3)*(a[15,1]+6*a[9,1]*a[4,1]+4*a[6,1]*a[7,1]+a[5,1]*a[11,1])/*
			 */-a[1,1]*(3*a[9,1]*a[3,1]+18*a[5,1]*a[3,1]*a[4,1]+12*a[2,1]*a[3,1]*a[7,1])/*
			 */-a[8,1]*(a[6,1]+6*a[2,1]*a[4,1]+a[11,1])
		f[,28]=ey2x5_d/*
			 */-(a[1,1]^2)*(a[15,1]+10*a[9,1]*a[4,1]+10*a[6,1]*a[7,1]+5*a[5,1]*a[11,1]/*
			 */+a[2,1]*a[14,1])/*
			 */-a[3,1]*(a[9,1]+10*a[5,1]*a[4,1]+10*a[2,1]*a[7,1]+a[14,1])
		f[,29]=eyx6_d/*
			 */-a[1,1]*(a[15,1]+15*a[9,1]*a[4,1]+20*a[6,1]*a[7,1]+15*a[5,1]*a[11,1]/*
			 */+6*a[2,1]*a[14,1])
	}

	f=f'
	return (f)
}


// subroutine to define the f vector for the optimal weighting matrix.
function optw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,x3_d,y3_d,y4x_d,/*
			*/y3x2_d,y2x3_d,yx4_d,y4_d,x4_d,y5x_d,y4x2_d,y3x3_d,y2x4_d,yx5_d,/*
            */y5_d,x5_d,y6x_d,y5x2_d,y4x3_d,y3x4_d,y2x5_d,yx6_d,/*
            */ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
            */ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
            */ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
            */ey3x4_d,ey2x5_d,eyx6_d,/*
            */inezz,zy_d,zx_d,eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,ey4_z,/*
            */ey3x_z,ey2x2_z,eyx3_z,ex4_z,ey5_z,ey4x_z,ey3x2_z,ey2x3_z,eyx4_z,ex5_z,/*
            */ey6_z,ey5x_z,ey4x2_z,ey3x3_z,ey2x4_z,eyx5_z,ex6_z,n,neq,estim)
{
	f=J(n,neq,0)

	f[,1]=y2_d:-ey2_d
	f[,2]=yx_d:-eyx_d
	f[,3]=x2_d:-ex2_d
	f[,4]=y2x_d:-ey2x_d
	f[,5]=yx2_d:-eyx2_d

	if (estim>1)
	{
		f[,6]=y3x_d:-ey3x_d
		f[,7]=y2x2_d:-ey2x2_d
		f[,8]=yx3_d:-eyx3_d
	}

	if (estim>2)
	{
		f[,9]=x3_d:-ex3_d
		f[,10]=y3_d:-ey3_d
		f[,11]=y4x_d:-ey4x_d
		f[,12]=y3x2_d:-ey3x2_d
		f[,13]=y2x3_d:-ey2x3_d
		f[,14]=yx4_d:-eyx4_d
	}

	if (estim>3)
	{
		f[,15]=y4_d:-ey4_d
		f[,16]=x4_d:-ex4_d
		f[,17]=y5x_d:-ey5x_d
		f[,18]=y4x2_d:-ey4x2_d
		f[,19]=y3x3_d:-ey3x3_d
		f[,20]=y2x4_d:-ey2x4_d
		f[,21]=yx5_d:-eyx5_d
	}

	if (estim>4)
	{
		f[,22]=y5_d:-ey5_d
		f[,23]=x5_d:-ex5_d
		f[,24]=y6x_d:-ey6x_d
		f[,25]=y5x2_d:-ey5x2_d
		f[,26]=y4x3_d:-ey4x3_d
		f[,27]=y3x4_d:-ey3x4_d
		f[,28]=y2x5_d:-ey2x5_d
		f[,29]=yx6_d:-eyx6_d
	}

	// this part makes the standard error adjustment.

	ei=J(n,neq,0)

	ei[,4]=(-2*eyx_z'*inezz*(zy_d')-ey2_z'*inezz*(zx_d'))'
	ei[,5]=(-ex2_z'*inezz*zy_d'-2*eyx_z'*inezz*zx_d')'

	if (estim>1)
	{
		ei[,6]=(-3*ey2x_z'*inezz*zy_d'-ey3_z'*inezz*zx_d')'
		ei[,7]=(-2*eyx2_z'*inezz*zy_d'-2*ey2x_z'*inezz*zx_d')'
		ei[,8]=(-ex3_z'*inezz*zy_d'-3*eyx2_z'*inezz*zx_d')'
	}

	if (estim>2)
	{
		ei[,9]=(-3*ex2_z'*inezz*zx_d')'
		ei[,10]=(-3*ey2_z'*inezz*zy_d')'
		ei[,11]=(-4*ey3x_z'*inezz*zy_d'-ey4_z'*inezz*zx_d')'
		ei[,12]=(-3*ey2x2_z'*inezz*zy_d'-2*ey3x_z'*inezz*zx_d')'
		ei[,13]=(-2*eyx3_z'*inezz*zy_d'-3*ey2x2_z'*inezz*zx_d')'
		ei[,14]=(-ex4_z'*inezz*zy_d'-4*eyx3_z'*inezz*zx_d')'
	}

	if (estim>3)
	{
		ei[,15]=(-4*ey3_z'*inezz*zy_d')'
		ei[,16]=(-4*ex3_z'*inezz*zx_d')'
		ei[,17]=(-5*ey4x_z'*inezz*zy_d'-ey5_z'*inezz*zx_d')'
		ei[,18]=(-4*ey3x2_z'*inezz*zy_d'-2*ey4x_z'*inezz*zx_d')'
		ei[,19]=(-3*ey2x3_z'*inezz*zy_d'-3*ey3x2_z'*inezz*zx_d')'
		ei[,20]=(-2*eyx4_z'*inezz*zy_d'-4*ey2x3_z'*inezz*zx_d')'
		ei[,21]=(-ex5_z'*inezz*zy_d'-5*eyx4_z'*inezz*zx_d')'
	}

	if (estim>4)
	{
		ei[,22]=(-5*ey4_z'*inezz*zy_d')'
		ei[,23]=(-5*ex4_z'*inezz*zx_d')'
		ei[,24]=(-6*ey5x_z'*inezz*zy_d'-ey6_z'*inezz*zx_d')'
		ei[,25]=(-5*ey4x2_z'*inezz*zy_d'-2*ey5x_z'*inezz*zx_d')'
		ei[,26]=(-4*ey3x3_z'*inezz*zy_d'-3*ey4x2_z'*inezz*zx_d')'
		ei[,27]=(-3*ey2x4_z'*inezz*zy_d'-4*ey3x3_z'*inezz*zx_d')'
		ei[,28]=(-2*eyx5_z'*inezz*zy_d'-5*ey2x4_z'*inezz*zx_d')'
		ei[,29]=(-ex6_z'*inezz*zy_d'-6*eyx5_z'*inezz*zx_d')'
	}

	f=f+ei

	return (f)
}


// subroutine to compute gradient matrix.
function grad1(a,neq)
{
	dvdc1=J(neq,1,0)
	dvdc2=J(neq,1,0)
	dvdc3=J(neq,1,0)
	dvdc4=J(neq,1,0)
	dvdc5=J(neq,1,0)

	dvdc1[1,1]=-2*a[1,1]*a[2,1]
	dvdc1[2,1]=-a[2,1]
	dvdc1[3,1]=0
	dvdc1[4,1]=-2*a[1,1]*a[5,1]
	dvdc1[5,1]=-a[5,1]

	dvdc2[1,1]=-(a[1,1]^2)
	dvdc2[2,1]=-a[1,1]
	dvdc2[3,1]=-1
	dvdc2[4,1]=0
	dvdc2[5,1]=0

	dvdc3[1,1]=-1
	dvdc3[2,1]=0
	dvdc3[3,1]=0
	dvdc3[4,1]=0
	dvdc3[5,1]=0

	dvdc4[1,1]=0
	dvdc4[2,1]=0
	dvdc4[3,1]=-1
	dvdc4[4,1]=0
	dvdc4[5,1]=0

	dvdc5[1,1]=0
	dvdc5[2,1]=0
	dvdc5[3,1]=0
	dvdc5[4,1]=-(a[1,1]^2)
	dvdc5[5,1]=-a[1,1]

	dvdc=(dvdc1,dvdc2,dvdc3,dvdc4,dvdc5)
	return(dvdc)
}

// subroutine to compute gradient matrix.
function grad2(a,neq)
{
	dvdc1=J(neq,1,0)
	dvdc2=J(neq,1,0)
	dvdc3=J(neq,1,0)
	dvdc4=J(neq,1,0)
	dvdc5=J(neq,1,0)
	dvdc6=J(neq,1,0)

	dvdc1[1,1]=-2*a[1,1]*a[2,1]
	dvdc1[2,1]=-a[2,1]
	dvdc1[3,1]=0
	dvdc1[4,1]=-2*a[1,1]*a[5,1]
	dvdc1[5,1]=-a[5,1]
	dvdc1[6,1]=-3*(a[1,1]^2)*a[6,1]-3*a[2,1]*a[3,1]
	dvdc1[7,1]=-2*a[1,1]*(a[6,1]+a[2,1]*a[4,1])
	dvdc1[8,1]=-a[6,1]-3*a[2,1]*a[4,1]

	dvdc2[1,1]=-(a[1,1]^2)
	dvdc2[2,1]=-a[1,1]
	dvdc2[3,1]=-1
	dvdc2[4,1]=0
	dvdc2[5,1]=0
	dvdc2[6,1]=-3*a[1,1]*a[3,1]
	dvdc2[7,1]=-(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc2[8,1]=-3*a[1,1]*a[4,1]

	dvdc3[1,1]=-1
	dvdc3[2,1]=0
	dvdc3[3,1]=0
	dvdc3[4,1]=0
	dvdc3[5,1]=0
	dvdc3[6,1]=-3*a[1,1]*a[2,1]
	dvdc3[7,1]=-(a[2,1]+a[4,1])
	dvdc3[8,1]=0

	dvdc4[1,1]=0
	dvdc4[2,1]=0
	dvdc4[3,1]=-1
	dvdc4[4,1]=0
	dvdc4[5,1]=0
	dvdc4[6,1]=0
	dvdc4[7,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc4[8,1]=-3*a[1,1]*a[2,1]

	dvdc5[1,1]=0
	dvdc5[2,1]=0
	dvdc5[3,1]=0
	dvdc5[4,1]=-(a[1,1]^2)
	dvdc5[5,1]=-a[1,1]
	dvdc5[6,1]=0
	dvdc5[7,1]=0
	dvdc5[8,1]=0

	dvdc6[1,1]=0
	dvdc6[2,1]=0
	dvdc6[3,1]=0
	dvdc6[4,1]=0
	dvdc6[5,1]=0
	dvdc6[6,1]=-(a[1,1]^3)
	dvdc6[7,1]=-(a[1,1]^2)
	dvdc6[8,1]=-a[1,1]

	dvdc=(dvdc1,dvdc2,dvdc3,dvdc4,dvdc5,dvdc6)
	return(dvdc)
}

// subroutine to compute gradient matrix.
function grad3(a,neq)
{
	dvdc1=J(neq,1,0)
	dvdc2=J(neq,1,0)
	dvdc3=J(neq,1,0)
	dvdc4=J(neq,1,0)
	dvdc5=J(neq,1,0)
	dvdc6=J(neq,1,0)
	dvdc7=J(neq,1,0)
	dvdc8=J(neq,1,0)
	dvdc9=J(neq,1,0)

	dvdc1[1,1]=-2*a[1,1]*a[2,1]
	dvdc1[2,1]=-a[2,1]
	dvdc1[3,1]=0
	dvdc1[4,1]=-2*a[1,1]*a[5,1]
	dvdc1[5,1]=-a[5,1]
	dvdc1[6,1]=-3*(a[1,1]^2)*a[6,1]-3*a[2,1]*a[3,1]
	dvdc1[7,1]=-2*a[1,1]*(a[6,1]+a[2,1]*a[4,1])
	dvdc1[8,1]=-a[6,1]-3*a[2,1]*a[4,1]
	dvdc1[9,1]=0
	dvdc1[10,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc1[11,1]=-4*(a[1,1]^3)*a[9,1]-12*a[1,1]*a[5,1]*a[3,1]/*
			  */-4*a[2,1]*a[8,1]
	dvdc1[12,1]=-3*(a[1,1]^2)*(a[9,1]+a[5,1]*a[4,1])/*
			  */-3*a[5,1]*a[3,1]
	dvdc1[13,1]=-2*a[1,1]*(a[9,1]+3*a[5,1]*a[4,1]+a[2,1]*a[7,1])
	dvdc1[14,1]=-(a[9,1]+6*a[5,1]*a[4,1]+4*a[2,1]*a[7,1])

	dvdc2[1,1]=-(a[1,1]^2)
	dvdc2[2,1]=-a[1,1]
	dvdc2[3,1]=-1
	dvdc2[4,1]=0
	dvdc2[5,1]=0
	dvdc2[6,1]=-3*a[1,1]*a[3,1]
	dvdc2[7,1]=-(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc2[8,1]=-3*a[1,1]*a[4,1]
	dvdc2[9,1]=0
	dvdc2[10,1]=0
	dvdc2[11,1]=-4*a[1,1]*a[8,1]
	dvdc2[12,1]=-a[8,1]
	dvdc2[13,1]=-(a[1,1]^2)*a[7,1]
	dvdc2[14,1]=-4*a[1,1]*a[7,1]

	dvdc3[1,1]=-1
	dvdc3[2,1]=0
	dvdc3[3,1]=0
	dvdc3[4,1]=0
	dvdc3[5,1]=0
	dvdc3[6,1]=-3*a[1,1]*a[2,1]
	dvdc3[7,1]=-(a[2,1]+a[4,1])
	dvdc3[8,1]=0
	dvdc3[9,1]=0
	dvdc3[10,1]=0
	dvdc3[11,1]=-6*(a[1,1]^2)*a[5,1]
	dvdc3[12,1]=-3*a[1,1]*a[5,1]
	dvdc3[13,1]=-(a[5,1]+a[7,1])
	dvdc3[14,1]=0

	dvdc4[1,1]=0
	dvdc4[2,1]=0
	dvdc4[3,1]=-1
	dvdc4[4,1]=0
	dvdc4[5,1]=0
	dvdc4[6,1]=0
	dvdc4[7,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc4[8,1]=-3*a[1,1]*a[2,1]
	dvdc4[9,1]=0
	dvdc4[10,1]=0
	dvdc4[11,1]=0
	dvdc4[12,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc4[13,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc4[14,1]=-6*a[1,1]*a[5,1]

	dvdc5[1,1]=0
	dvdc5[2,1]=0
	dvdc5[3,1]=0
	dvdc5[4,1]=-(a[1,1]^2)
	dvdc5[5,1]=-a[1,1]
	dvdc5[6,1]=0
	dvdc5[7,1]=0
	dvdc5[8,1]=0
	dvdc5[9,1]=-1
	dvdc5[10,1]=-(a[1,1]^3)
	dvdc5[11,1]=-6*(a[1,1]^2)*a[3,1]
	dvdc5[12,1]=-(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc5[13,1]=-3*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc5[14,1]=-6*a[1,1]*a[4,1]

	dvdc6[1,1]=0
	dvdc6[2,1]=0
	dvdc6[3,1]=0
	dvdc6[4,1]=0
	dvdc6[5,1]=0
	dvdc6[6,1]=-(a[1,1]^3)
	dvdc6[7,1]=-(a[1,1]^2)
	dvdc6[8,1]=-a[1,1]
	dvdc6[9,1]=0
	dvdc6[10,1]=0
	dvdc6[11,1]=0
	dvdc6[12,1]=0
	dvdc6[13,1]=0
	dvdc6[14,1]=0

	dvdc7[1,1]=0
	dvdc7[2,1]=0
	dvdc7[3,1]=0
	dvdc7[4,1]=0
	dvdc7[5,1]=0
	dvdc7[6,1]=0
	dvdc7[7,1]=0
	dvdc7[8,1]=0
	dvdc7[9,1]=-1
	dvdc7[10,1]=0
	dvdc7[11,1]=0
	dvdc7[12,1]=0
	dvdc7[13,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc7[14,1]=-4*a[1,1]*a[2,1]

	dvdc8[1,1]=0
	dvdc8[2,1]=0
	dvdc8[3,1]=0
	dvdc8[4,1]=0
	dvdc8[5,1]=0
	dvdc8[6,1]=0
	dvdc8[7,1]=0
	dvdc8[8,1]=0
	dvdc8[9,1]=0
	dvdc8[10,1]=-1
	dvdc8[11,1]=-4*a[1,1]*a[2,1]
	dvdc8[12,1]=-(a[2,1]+a[4,1])
	dvdc8[13,1]=0
	dvdc8[14,1]=0

	dvdc9[1,1]=0
	dvdc9[2,1]=0
	dvdc9[3,1]=0
	dvdc9[4,1]=0
	dvdc9[5,1]=0
	dvdc9[6,1]=0
	dvdc9[7,1]=0
	dvdc9[8,1]=0
	dvdc9[9,1]=0
	dvdc9[10,1]=0
	dvdc9[11,1]=-(a[1,1]^4)
	dvdc9[12,1]=-(a[1,1]^3)
	dvdc9[13,1]=-(a[1,1]^2)
	dvdc9[14,1]=-a[1,1]

	dvdc=(dvdc1,dvdc2,dvdc3,dvdc4,dvdc5,dvdc6,dvdc7,dvdc8,dvdc9)
	return(dvdc)
}



// subroutine to compute gradient matrix.
function grad4(a,neq)
{
	dvdc1=J(neq,1,0)
	dvdc2=J(neq,1,0)
	dvdc3=J(neq,1,0)
	dvdc4=J(neq,1,0)
	dvdc5=J(neq,1,0)
	dvdc6=J(neq,1,0)
	dvdc7=J(neq,1,0)
	dvdc8=J(neq,1,0)
	dvdc9=J(neq,1,0)
	dvdc10=J(neq,1,0)
	dvdc11=J(neq,1,0)
	dvdc12=J(neq,1,0)

	dvdc1[1,1]=-2*a[1,1]*a[2,1]
	dvdc1[2,1]=-a[2,1]
	dvdc1[3,1]=0
	dvdc1[4,1]=-2*a[1,1]*a[5,1]
	dvdc1[5,1]=-a[5,1]
	dvdc1[6,1]=-3*(a[1,1]^2)*a[6,1]-3*a[2,1]*a[3,1]
	dvdc1[7,1]=-2*a[1,1]*(a[6,1]+a[2,1]*a[4,1])
	dvdc1[8,1]=-a[6,1]-3*a[2,1]*a[4,1]
	dvdc1[9,1]=0
	dvdc1[10,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc1[11,1]=-4*(a[1,1]^3)*a[9,1]-12*a[1,1]*a[5,1]*a[3,1]/*
			  */-4*a[2,1]*a[8,1]
	dvdc1[12,1]=-3*(a[1,1]^2)*(a[9,1]+a[5,1]*a[4,1])/*
			  */-3*a[5,1]*a[3,1]
	dvdc1[13,1]=-2*a[1,1]*(a[9,1]+3*a[5,1]*a[4,1]+a[2,1]*a[7,1])
	dvdc1[14,1]=-(a[9,1]+6*a[5,1]*a[4,1]+4*a[2,1]*a[7,1])
	dvdc1[15,1]=-4*(a[1,1]^3)*a[6,1]-12*a[1,1]*a[2,1]*a[3,1]
	dvdc1[16,1]=0
	dvdc1[17,1]=-5*(a[1,1]^4)*a[12,1]-30*(a[1,1]^2)*a[6,1]*a[3,1]/*
			  */-20*a[1,1]*a[5,1]*a[8,1]-5*a[2,1]*a[10,1]
	dvdc1[18,1]=-4*(a[1,1]^3)*(a[12,1]+a[6,1]*a[4,1])/*
			  */-12*a[1,1]*(a[6,1]*a[3,1]+a[2,1]*a[3,1]*a[4,1])/*
			  */-4*a[5,1]*a[8,1]
	dvdc1[19,1]=-3*(a[1,1]^2)*(a[12,1]+3*a[6,1]*a[4,1]+a[5,1]*a[7,1])/*
			  */-(3*a[6,1]*a[3,1]+9*a[2,1]*a[3,1]*a[4,1])
	dvdc1[20,1]=-2*a[1,1]*(a[12,1]+6*a[6,1]*a[4,1]+4*a[5,1]*a[7,1]/*
			  */+a[2,1]*a[11,1])
	dvdc1[21,1]=-(a[12,1]+10*a[6,1]*a[4,1]+10*a[5,1]*a[7,1]/*
			  */+5*a[2,1]*a[11,1])

	dvdc2[1,1]=-(a[1,1]^2)
	dvdc2[2,1]=-a[1,1]
	dvdc2[3,1]=-1
	dvdc2[4,1]=0
	dvdc2[5,1]=0
	dvdc2[6,1]=-3*a[1,1]*a[3,1]
	dvdc2[7,1]=-(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc2[8,1]=-3*a[1,1]*a[4,1]
	dvdc2[9,1]=0
	dvdc2[10,1]=0
	dvdc2[11,1]=-4*a[1,1]*a[8,1]
	dvdc2[12,1]=-a[8,1]
	dvdc2[13,1]=-(a[1,1]^2)*a[7,1]
	dvdc2[14,1]=-4*a[1,1]*a[7,1]
	dvdc2[15,1]=-6*(a[1,1]^2)*a[3,1]
	dvdc2[16,1]=-6*a[4,1]
	dvdc2[17,1]=-5*a[1,1]*a[10,1]
	dvdc2[18,1]=-6*(a[1,1]^2)*a[3,1]*a[4,1]-a[10,1]
	dvdc2[19,1]=-9*a[1,1]*a[3,1]*a[4,1]
	dvdc2[20,1]=-(a[1,1]^2)*a[11,1]-6*a[3,1]*a[4,1]
	dvdc2[21,1]=-5*a[1,1]*a[11,1]

	dvdc3[1,1]=-1
	dvdc3[2,1]=0
	dvdc3[3,1]=0
	dvdc3[4,1]=0
	dvdc3[5,1]=0
	dvdc3[6,1]=-3*a[1,1]*a[2,1]
	dvdc3[7,1]=-(a[2,1]+a[4,1])
	dvdc3[8,1]=0
	dvdc3[9,1]=0
	dvdc3[10,1]=0
	dvdc3[11,1]=-6*(a[1,1]^2)*a[5,1]
	dvdc3[12,1]=-3*a[1,1]*a[5,1]
	dvdc3[13,1]=-(a[5,1]+a[7,1])
	dvdc3[14,1]=0
	dvdc3[15,1]=-6*(a[1,1]^2)*a[2,1]
	dvdc3[16,1]=0
	dvdc3[17,1]=-10*(a[1,1]^3)*a[6,1]
	dvdc3[18,1]=-6*(a[1,1]^2)*(a[6,1]+a[2,1])
	dvdc3[19,1]=-a[1,1]*(3*a[6,1]+9*a[2,1]*a[4,1])
	dvdc3[20,1]=-(a[6,1]+6*a[2,1]*a[4,1]+a[11,1])
	dvdc3[21,1]=0

	dvdc4[1,1]=0
	dvdc4[2,1]=0
	dvdc4[3,1]=-1
	dvdc4[4,1]=0
	dvdc4[5,1]=0
	dvdc4[6,1]=0
	dvdc4[7,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc4[8,1]=-3*a[1,1]*a[2,1]
	dvdc4[9,1]=0
	dvdc4[10,1]=0
	dvdc4[11,1]=0
	dvdc4[12,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc4[13,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc4[14,1]=-6*a[1,1]*a[5,1]
	dvdc4[15,1]=0
	dvdc4[16,1]=-6*a[2,1]
	dvdc4[17,1]=0
	dvdc4[18,1]=-(a[1,1]^4)*a[6,1]-6*(a[1,1]^2)*a[2,1]*a[3,1]-a[10,1]
	dvdc4[19,1]=-3*(a[1,1]^3)*a[6,1]-9*a[1,1]*a[2,1]*a[3,1]
	dvdc4[20,1]=-6*(a[1,1]^2)*a[6,1]-6*a[2,1]*a[3,1]
	dvdc4[21,1]=-10*a[1,1]*a[6,1]

	dvdc5[1,1]=0
	dvdc5[2,1]=0
	dvdc5[3,1]=0
	dvdc5[4,1]=-(a[1,1]^2)
	dvdc5[5,1]=-a[1,1]
	dvdc5[6,1]=0
	dvdc5[7,1]=0
	dvdc5[8,1]=0
	dvdc5[9,1]=-1
	dvdc5[10,1]=-(a[1,1]^3)
	dvdc5[11,1]=-6*(a[1,1]^2)*a[3,1]
	dvdc5[12,1]=-(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc5[13,1]=-3*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc5[14,1]=-6*a[1,1]*a[4,1]
	dvdc5[15,1]=0
	dvdc5[16,1]=0
	dvdc5[17,1]=-10*(a[1,1]^2)*a[8,1]
	dvdc5[18,1]=-4*a[1,1]*a[8,1]
	dvdc5[19,1]=-(a[1,1]^3)*a[7,1]-a[8,1]
	dvdc5[20,1]=-4*(a[1,1]^2)*a[7,1]
	dvdc5[21,1]=-10*a[1,1]*a[7,1]

	dvdc6[1,1]=0
	dvdc6[2,1]=0
	dvdc6[3,1]=0
	dvdc6[4,1]=0
	dvdc6[5,1]=0
	dvdc6[6,1]=-(a[1,1]^3)
	dvdc6[7,1]=-(a[1,1]^2)
	dvdc6[8,1]=-a[1,1]
	dvdc6[9,1]=0
	dvdc6[10,1]=0
	dvdc6[11,1]=0
	dvdc6[12,1]=0
	dvdc6[13,1]=0
	dvdc6[14,1]=0
	dvdc6[15,1]=-(a[1,1]^4)
	dvdc6[16,1]=-1
	dvdc6[17,1]=-10*(a[1,1]^3)*a[3,1]
	dvdc6[18,1]=-(a[1,1]^4)*a[4,1]-6*(a[1,1]^2)*a[3,1]
	dvdc6[19,1]=-3*(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc6[20,1]=-6*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc6[21,1]=-10*a[1,1]*a[4,1]

	dvdc7[1,1]=0
	dvdc7[2,1]=0
	dvdc7[3,1]=0
	dvdc7[4,1]=0
	dvdc7[5,1]=0
	dvdc7[6,1]=0
	dvdc7[7,1]=0
	dvdc7[8,1]=0
	dvdc7[9,1]=-1
	dvdc7[10,1]=0
	dvdc7[11,1]=0
	dvdc7[12,1]=0
	dvdc7[13,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc7[14,1]=-4*a[1,1]*a[2,1]
	dvdc7[15,1]=0
	dvdc7[16,1]=0
	dvdc7[17,1]=0
	dvdc7[18,1]=0
	dvdc7[19,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc7[20,1]=-4*(a[1,1]^2)*a[5,1]
	dvdc7[21,1]=-10*a[1,1]*a[5,1]

	dvdc8[1,1]=0
	dvdc8[2,1]=0
	dvdc8[3,1]=0
	dvdc8[4,1]=0
	dvdc8[5,1]=0
	dvdc8[6,1]=0
	dvdc8[7,1]=0
	dvdc8[8,1]=0
	dvdc8[9,1]=0
	dvdc8[10,1]=-1
	dvdc8[11,1]=-4*a[1,1]*a[2,1]
	dvdc8[12,1]=-(a[2,1]+a[4,1])
	dvdc8[13,1]=0
	dvdc8[14,1]=0
	dvdc8[15,1]=0
	dvdc8[16,1]=0
	dvdc8[17,1]=-10*(a[1,1]^2)*a[5,1]
	dvdc8[18,1]=-4*a[1,1]*a[5,1]
	dvdc8[19,1]=-(a[5,1]+a[7,1])
	dvdc8[20,1]=0
	dvdc8[21,1]=0

	dvdc9[1,1]=0
	dvdc9[2,1]=0
	dvdc9[3,1]=0
	dvdc9[4,1]=0
	dvdc9[5,1]=0
	dvdc9[6,1]=0
	dvdc9[7,1]=0
	dvdc9[8,1]=0
	dvdc9[9,1]=0
	dvdc9[10,1]=0
	dvdc9[11,1]=-(a[1,1]^4)
	dvdc9[12,1]=-(a[1,1]^3)
	dvdc9[13,1]=-(a[1,1]^2)
	dvdc9[14,1]=-a[1,1]
	dvdc9[15,1]=0
	dvdc9[16,1]=0
	dvdc9[17,1]=0
	dvdc9[18,1]=0
	dvdc9[19,1]=0
	dvdc9[20,1]=0
	dvdc9[21,1]=0

	dvdc10[1,1]=0
	dvdc10[2,1]=0
	dvdc10[3,1]=0
	dvdc10[4,1]=0
	dvdc10[5,1]=0
	dvdc10[6,1]=0
	dvdc10[7,1]=0
	dvdc10[8,1]=0
	dvdc10[9,1]=0
	dvdc10[10,1]=0
	dvdc10[11,1]=0
	dvdc10[12,1]=0
	dvdc10[13,1]=0
	dvdc10[14,1]=0
	dvdc10[15,1]=-1
	dvdc10[16,1]=0
	dvdc10[17,1]=-5*a[1,1]*a[2,1]
	dvdc10[18,1]=-(a[2,1]+a[4,1])
	dvdc10[19,1]=0
	dvdc10[20,1]=0
	dvdc10[21,1]=0

	dvdc11[1,1]=0
	dvdc11[2,1]=0
	dvdc11[3,1]=0
	dvdc11[4,1]=0
	dvdc11[5,1]=0
	dvdc11[6,1]=0
	dvdc11[7,1]=0
	dvdc11[8,1]=0
	dvdc11[9,1]=0
	dvdc11[10,1]=0
	dvdc11[11,1]=0
	dvdc11[12,1]=0
	dvdc11[13,1]=0
	dvdc11[14,1]=0
	dvdc11[15,1]=0
	dvdc11[16,1]=-1
	dvdc11[17,1]=0
	dvdc11[18,1]=0
	dvdc11[19,1]=0
	dvdc11[20,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc11[21,1]=-5*a[1,1]*a[2,1]

	dvdc12[1,1]=0
	dvdc12[2,1]=0
	dvdc12[3,1]=0
	dvdc12[4,1]=0
	dvdc12[5,1]=0
	dvdc12[6,1]=0
	dvdc12[7,1]=0
	dvdc12[8,1]=0
	dvdc12[9,1]=0
	dvdc12[10,1]=0
	dvdc12[11,1]=0
	dvdc12[12,1]=0
	dvdc12[13,1]=0
	dvdc12[14,1]=0
	dvdc12[15,1]=0
	dvdc12[16,1]=0
	dvdc12[17,1]=-(a[1,1]^5)
	dvdc12[18,1]=-(a[1,1]^4)
	dvdc12[19,1]=-(a[1,1]^3)
	dvdc12[20,1]=-(a[1,1]^2)
	dvdc12[21,1]=-a[1,1]

	dvdc=(dvdc1,dvdc2,dvdc3,dvdc4,dvdc5,dvdc6,dvdc7,dvdc8,dvdc9,dvdc10,dvdc11,dvdc12)
	return(dvdc)
}



// subroutine to compute gradient matrix.
function grad5(a,neq)
{
	dvdc1=J(neq,1,0)
	dvdc2=J(neq,1,0)
	dvdc3=J(neq,1,0)
	dvdc4=J(neq,1,0)
	dvdc5=J(neq,1,0)
	dvdc6=J(neq,1,0)
	dvdc7=J(neq,1,0)
	dvdc8=J(neq,1,0)
	dvdc9=J(neq,1,0)
	dvdc10=J(neq,1,0)
	dvdc11=J(neq,1,0)
	dvdc12=J(neq,1,0)
	dvdc13=J(neq,1,0)
	dvdc14=J(neq,1,0)
	dvdc15=J(neq,1,0)

	dvdc1[1,1]=-2*a[1,1]*a[2,1]
	dvdc1[2,1]=-a[2,1]
	dvdc1[3,1]=0
	dvdc1[4,1]=-2*a[1,1]*a[5,1]
	dvdc1[5,1]=-a[5,1]
	dvdc1[6,1]=-3*(a[1,1]^2)*a[6,1]-3*a[2,1]*a[3,1]
	dvdc1[7,1]=-2*a[1,1]*(a[6,1]+a[2,1]*a[4,1])
	dvdc1[8,1]=-a[6,1]-3*a[2,1]*a[4,1]
	dvdc1[9,1]=0
	dvdc1[10,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc1[11,1]=-4*(a[1,1]^3)*a[9,1]-12*a[1,1]*a[5,1]*a[3,1]/*
			  */-4*a[2,1]*a[8,1]
	dvdc1[12,1]=-3*(a[1,1]^2)*(a[9,1]+a[5,1]*a[4,1])/*
			  */-3*a[5,1]*a[3,1]
	dvdc1[13,1]=-2*a[1,1]*(a[9,1]+3*a[5,1]*a[4,1]+a[2,1]*a[7,1])
	dvdc1[14,1]=-(a[9,1]+6*a[5,1]*a[4,1]+4*a[2,1]*a[7,1])
	dvdc1[15,1]=-4*(a[1,1]^3)*a[6,1]-12*a[1,1]*a[2,1]*a[3,1]
	dvdc1[16,1]=0
	dvdc1[17,1]=-5*(a[1,1]^4)*a[12,1]-30*(a[1,1]^2)*a[6,1]*a[3,1]/*
			  */-20*a[1,1]*a[5,1]*a[8,1]-5*a[2,1]*a[10,1]
	dvdc1[18,1]=-4*(a[1,1]^3)*(a[12,1]+a[6,1]*a[4,1])/*
			  */-12*a[1,1]*(a[6,1]*a[3,1]+a[2,1]*a[3,1]*a[4,1])/*
			  */-4*a[5,1]*a[8,1]
	dvdc1[19,1]=-3*(a[1,1]^2)*(a[12,1]+3*a[6,1]*a[4,1]+a[5,1]*a[7,1])/*
			  */-(3*a[6,1]*a[3,1]+9*a[2,1]*a[3,1]*a[4,1])
	dvdc1[20,1]=-2*a[1,1]*(a[12,1]+6*a[6,1]*a[4,1]+4*a[5,1]*a[7,1]/*
			  */+a[2,1]*a[11,1])
	dvdc1[21,1]=-(a[12,1]+10*a[6,1]*a[4,1]+10*a[5,1]*a[7,1]/*
			  */+5*a[2,1]*a[11,1])
	dvdc1[22,1]=-5*(a[1,1]^4)*a[9,1]-30*(a[1,1]^2)*a[5,1]*a[3,1]/*
			  */-20*a[1,1]*a[2,1]*a[8,1]
	dvdc1[23,1]=0
	dvdc1[24,1]=-6*(a[1,1]^5)*a[15,1]-60*(a[1,1]^3)*a[9,1]*a[3,1]/*
			  */-60*(a[1,1]^2)*a[6,1]*a[8,1]-30*a[1,1]*a[5,1]*a[10,1]/*
			  */-6*a[2,1]*a[13,1]
	dvdc1[25,1]=-5*(a[1,1]^4)*(a[15,1]+a[9,1]*a[4,1])/*
			  */-30*(a[1,1]^2)*(a[9,1]*a[3,1]+a[5,1]*a[3,1]*a[4,1])/*
			  */-20*a[1,1]*(a[6,1]*a[8,1]+a[2,1]*a[8,1]*a[4,1])/*
			  */-5*a[5,1]*a[10,1]
	dvdc1[26,1]=-4*(a[1,1]^3)*(a[15,1]+3*a[9,1]*a[4,1]+a[6,1]*a[7,1])/*
			  */-12*a[1,1]*(a[9,1]*a[3,1]+3*a[5,1]*a[3,1]*a[4,1]+a[2,1]*a[3,1]*a[7,1])/*
			  */-4*(a[6,1]*a[8,1]+3*a[2,1]*a[8,1]*a[4,1])
	dvdc1[27,1]=-3*(a[1,1]^2)*(a[15,1]+6*a[9,1]*a[4,1]+4*a[6,1]*a[7,1]+a[5,1]*a[11,1])/*
			  */-(3*a[9,1]*a[3,1]+18*a[5,1]*a[3,1]*a[4,1]+12*a[2,1]*a[3,1]*a[7,1])
	dvdc1[28,1]=-2*a[1,1]*(a[15,1]+10*a[9,1]*a[4,1]+10*a[6,1]*a[7,1]+5*a[5,1]*a[11,1]/*
			  */+a[2,1]*a[14,1])
	dvdc1[29,1]=-(a[15,1]+15*a[9,1]*a[4,1]+20*a[6,1]*a[7,1]/*
			  */+15*a[5,1]*a[11,1]+6*a[2,1]*a[14,1])

	dvdc2[1,1]=-(a[1,1]^2)
	dvdc2[2,1]=-a[1,1]
	dvdc2[3,1]=-1
	dvdc2[4,1]=0
	dvdc2[5,1]=0
	dvdc2[6,1]=-3*a[1,1]*a[3,1]
	dvdc2[7,1]=-(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc2[8,1]=-3*a[1,1]*a[4,1]
	dvdc2[9,1]=0
	dvdc2[10,1]=0
	dvdc2[11,1]=-4*a[1,1]*a[8,1]
	dvdc2[12,1]=-a[8,1]
	dvdc2[13,1]=-(a[1,1]^2)*a[7,1]
	dvdc2[14,1]=-4*a[1,1]*a[7,1]
	dvdc2[15,1]=-6*(a[1,1]^2)*a[3,1]
	dvdc2[16,1]=-6*a[4,1]
	dvdc2[17,1]=-5*a[1,1]*a[10,1]
	dvdc2[18,1]=-6*(a[1,1]^2)*a[3,1]*a[4,1]-a[10,1]
	dvdc2[19,1]=-9*a[1,1]*a[3,1]*a[4,1]
	dvdc2[20,1]=-(a[1,1]^2)*a[11,1]-6*a[3,1]*a[4,1]
	dvdc2[21,1]=-5*a[1,1]*a[11,1]
	dvdc2[22,1]=-10*(a[1,1]^2)*a[8,1]
	dvdc2[23,1]=-10*a[7,1]
	dvdc2[24,1]=-6*a[1,1]*a[13,1]
	dvdc2[25,1]=-10*(a[1,1]^2)*a[8,1]*a[4,1]-a[13,1]
	dvdc2[26,1]=-6*(a[1,1]^2)*a[3,1]*a[7,1]-12*a[1,1]*a[8,1]*a[4,1]
	dvdc2[27,1]=-12*a[1,1]*a[3,1]*a[7,1]-6*a[8,1]*a[4,1]
	dvdc2[28,1]=-(a[1,1]^2)*a[14,1]-10*a[3,1]*a[7,1]
	dvdc2[29,1]=-6*a[1,1]*a[14,1]

	dvdc3[1,1]=-1
	dvdc3[2,1]=0
	dvdc3[3,1]=0
	dvdc3[4,1]=0
	dvdc3[5,1]=0
	dvdc3[6,1]=-3*a[1,1]*a[2,1]
	dvdc3[7,1]=-(a[2,1]+a[4,1])
	dvdc3[8,1]=0
	dvdc3[9,1]=0
	dvdc3[10,1]=0
	dvdc3[11,1]=-6*(a[1,1]^2)*a[5,1]
	dvdc3[12,1]=-3*a[1,1]*a[5,1]
	dvdc3[13,1]=-(a[5,1]+a[7,1])
	dvdc3[14,1]=0
	dvdc3[15,1]=-6*(a[1,1]^2)*a[2,1]
	dvdc3[16,1]=0
	dvdc3[17,1]=-10*(a[1,1]^3)*a[6,1]
	dvdc3[18,1]=-6*(a[1,1]^2)*(a[6,1]+a[2,1])
	dvdc3[19,1]=-a[1,1]*(3*a[6,1]+9*a[2,1]*a[4,1])
	dvdc3[20,1]=-(a[6,1]+6*a[2,1]*a[4,1]+a[11,1])
	dvdc3[21,1]=0
	dvdc3[22,1]=-10*(a[1,1]^3)*a[5,1]
	dvdc3[23,1]=0
	dvdc3[24,1]=-15*(a[1,1]^4)*a[9,1]
	dvdc3[25,1]=-10*(a[1,1]^3)*(a[9,1]+a[5,1]*a[4,1])
	dvdc3[26,1]=-6*(a[1,1]^2)*(a[9,1]+3*a[5,1]*a[4,1]+a[2,1]*a[7,1])
	dvdc3[27,1]=-a[1,1]*(3*a[9,1]+18*a[5,1]*a[4,1]+12*a[2,1]*a[7,1])
	dvdc3[28,1]=-(a[9,1]+10*a[5,1]*a[4,1]+10*a[2,1]*a[7,1]+a[14,1])
	dvdc3[29,1]=0

	dvdc4[1,1]=0
	dvdc4[2,1]=0
	dvdc4[3,1]=-1
	dvdc4[4,1]=0
	dvdc4[5,1]=0
	dvdc4[6,1]=0
	dvdc4[7,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc4[8,1]=-3*a[1,1]*a[2,1]
	dvdc4[9,1]=0
	dvdc4[10,1]=0
	dvdc4[11,1]=0
	dvdc4[12,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc4[13,1]=-3*(a[1,1]^2)*a[5,1]
	dvdc4[14,1]=-6*a[1,1]*a[5,1]
	dvdc4[15,1]=0
	dvdc4[16,1]=-6*a[2,1]
	dvdc4[17,1]=0
	dvdc4[18,1]=-(a[1,1]^4)*a[6,1]-6*(a[1,1]^2)*a[2,1]*a[3,1]-a[10,1]
	dvdc4[19,1]=-3*(a[1,1]^3)*a[6,1]-9*a[1,1]*a[2,1]*a[3,1]
	dvdc4[20,1]=-6*(a[1,1]^2)*a[6,1]-6*a[2,1]*a[3,1]
	dvdc4[21,1]=-10*a[1,1]*a[6,1]
	dvdc4[22,1]=0
	dvdc4[23,1]=-10*a[5,1]
	dvdc4[24,1]=0
	dvdc4[25,1]=-(a[1,1]^5)*a[9,1]-10*(a[1,1]^3)*a[5,1]*a[3,1]/*
			  */-10*(a[1,1]^2)*a[2,1]*a[8,1]-a[13,1]
	dvdc4[26,1]=-3*(a[1,1]^4)*a[9,1]-18*(a[1,1]^2)*a[5,1]*a[3,1]/*
			  */-12*a[1,1]*a[2,1]*a[8,1]
	dvdc4[27,1]=-6*(a[1,1]^3)*a[9,1]-18*a[1,1]*a[5,1]*a[3,1]/*
			  */-6*a[8,1]*a[2,1]
	dvdc4[28,1]=-10*(a[1,1]^2)*a[9,1]-10*a[3,1]*a[5,1]
	dvdc4[29,1]=-15*a[1,1]*a[9,1]

	dvdc5[1,1]=0
	dvdc5[2,1]=0
	dvdc5[3,1]=0
	dvdc5[4,1]=-(a[1,1]^2)
	dvdc5[5,1]=-a[1,1]
	dvdc5[6,1]=0
	dvdc5[7,1]=0
	dvdc5[8,1]=0
	dvdc5[9,1]=-1
	dvdc5[10,1]=-(a[1,1]^3)
	dvdc5[11,1]=-6*(a[1,1]^2)*a[3,1]
	dvdc5[12,1]=-(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc5[13,1]=-3*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc5[14,1]=-6*a[1,1]*a[4,1]
	dvdc5[15,1]=0
	dvdc5[16,1]=0
	dvdc5[17,1]=-10*(a[1,1]^2)*a[8,1]
	dvdc5[18,1]=-4*a[1,1]*a[8,1]
	dvdc5[19,1]=-(a[1,1]^3)*a[7,1]-a[8,1]
	dvdc5[20,1]=-4*(a[1,1]^2)*a[7,1]
	dvdc5[21,1]=-10*a[1,1]*a[7,1]
	dvdc5[22,1]=-10*(a[1,1]^3)*a[3,1]
	dvdc5[23,1]=-10*a[4,1]
	dvdc5[24,1]=-15*(a[1,1]^2)*a[10,1]
	dvdc5[25,1]=-10*(a[1,1]^3)*a[3,1]*a[4,1]-5*a[1,1]*a[10,1]
	dvdc5[26,1]=-18*(a[1,1]^2)*a[3,1]*a[4,1]-a[10,1]
	dvdc5[27,1]=-(a[1,1]^3)*a[11,1]-18*a[1,1]*a[3,1]*a[4,1]
	dvdc5[28,1]=-5*(a[1,1]^2)*a[11,1]-10*a[3,1]*a[4,1]
	dvdc5[29,1]=-15*a[1,1]*a[11,1]

	dvdc6[1,1]=0
	dvdc6[2,1]=0
	dvdc6[3,1]=0
	dvdc6[4,1]=0
	dvdc6[5,1]=0
	dvdc6[6,1]=-(a[1,1]^3)
	dvdc6[7,1]=-(a[1,1]^2)
	dvdc6[8,1]=-a[1,1]
	dvdc6[9,1]=0
	dvdc6[10,1]=0
	dvdc6[11,1]=0
	dvdc6[12,1]=0
	dvdc6[13,1]=0
	dvdc6[14,1]=0
	dvdc6[15,1]=-(a[1,1]^4)
	dvdc6[16,1]=-1
	dvdc6[17,1]=-10*(a[1,1]^3)*a[3,1]
	dvdc6[18,1]=-(a[1,1]^4)*a[4,1]-6*(a[1,1]^2)*a[3,1]
	dvdc6[19,1]=-3*(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc6[20,1]=-6*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc6[21,1]=-10*a[1,1]*a[4,1]
	dvdc6[22,1]=0
	dvdc6[23,1]=0
	dvdc6[24,1]=-20*(a[1,1]^3)*a[8,1]
	dvdc6[25,1]=-10*(a[1,1]^2)*a[8,1]
	dvdc6[26,1]=-(a[1,1]^4)*a[7,1]-4*a[1,1]*a[8,1]
	dvdc6[27,1]=-4*(a[1,1]^3)*a[7,1]-a[8,1]
	dvdc6[28,1]=-10*(a[1,1]^2)*a[7,1]
	dvdc6[29,1]=-20*a[1,1]*a[7,1]

	dvdc7[1,1]=0
	dvdc7[2,1]=0
	dvdc7[3,1]=0
	dvdc7[4,1]=0
	dvdc7[5,1]=0
	dvdc7[6,1]=0
	dvdc7[7,1]=0
	dvdc7[8,1]=0
	dvdc7[9,1]=-1
	dvdc7[10,1]=0
	dvdc7[11,1]=0
	dvdc7[12,1]=0
	dvdc7[13,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc7[14,1]=-4*a[1,1]*a[2,1]
	dvdc7[15,1]=0
	dvdc7[16,1]=0
	dvdc7[17,1]=0
	dvdc7[18,1]=0
	dvdc7[19,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc7[20,1]=-4*(a[1,1]^2)*a[5,1]
	dvdc7[21,1]=-10*a[1,1]*a[5,1]
	dvdc7[22,1]=0
	dvdc7[23,1]=-10*a[2,1]
	dvdc7[24,1]=0
	dvdc7[25,1]=0
	dvdc7[26,1]=-(a[1,1]^4)*a[6,1]-6*(a[1,1]^2)*a[2,1]*a[3,1]-a[10,1]
	dvdc7[27,1]=-4*(a[1,1]^3)*a[6,1]-12*a[1,1]*a[2,1]*a[3,1]
	dvdc7[28,1]=-10*(a[1,1]^2)*a[6,1]-10*a[3,1]*a[2,1]
	dvdc7[29,1]=-20*a[1,1]*a[6,1]

	dvdc8[1,1]=0
	dvdc8[2,1]=0
	dvdc8[3,1]=0
	dvdc8[4,1]=0
	dvdc8[5,1]=0
	dvdc8[6,1]=0
	dvdc8[7,1]=0
	dvdc8[8,1]=0
	dvdc8[9,1]=0
	dvdc8[10,1]=-1
	dvdc8[11,1]=-4*a[1,1]*a[2,1]
	dvdc8[12,1]=-(a[2,1]+a[4,1])
	dvdc8[13,1]=0
	dvdc8[14,1]=0
	dvdc8[15,1]=0
	dvdc8[16,1]=0
	dvdc8[17,1]=-10*(a[1,1]^2)*a[5,1]
	dvdc8[18,1]=-4*a[1,1]*a[5,1]
	dvdc8[19,1]=-(a[5,1]+a[7,1])
	dvdc8[20,1]=0
	dvdc8[21,1]=0
	dvdc8[22,1]=-10*(a[1,1]^2)*a[2,1]
	dvdc8[23,1]=0
	dvdc8[24,1]=-20*(a[1,1]^3)*a[6,1]
	dvdc8[25,1]=-10*(a[1,1]^2)*(a[6,1]+a[2,1]*a[4,1])
	dvdc8[26,1]=-4*a[1,1]*(a[6,1]+3*a[2,1]*a[4,1])
	dvdc8[27,1]=-(a[6,1]+6*a[2,1]*a[4,1]+a[11,1])
	dvdc8[28,1]=0
	dvdc8[29,1]=0

	dvdc9[1,1]=0
	dvdc9[2,1]=0
	dvdc9[3,1]=0
	dvdc9[4,1]=0
	dvdc9[5,1]=0
	dvdc9[6,1]=0
	dvdc9[7,1]=0
	dvdc9[8,1]=0
	dvdc9[9,1]=0
	dvdc9[10,1]=0
	dvdc9[11,1]=-(a[1,1]^4)
	dvdc9[12,1]=-(a[1,1]^3)
	dvdc9[13,1]=-(a[1,1]^2)
	dvdc9[14,1]=-a[1,1]
	dvdc9[15,1]=0
	dvdc9[16,1]=0
	dvdc9[17,1]=0
	dvdc9[18,1]=0
	dvdc9[19,1]=0
	dvdc9[20,1]=0
	dvdc9[21,1]=0
	dvdc9[22,1]=-(a[1,1]^5)
	dvdc9[23,1]=-1
	dvdc9[24,1]=-15*(a[1,1]^4)*a[3,1]
	dvdc9[25,1]=-(a[1,1]^5)*a[4,1]-10*(a[1,1]^3)*a[3,1]
	dvdc9[26,1]=-3*(a[1,1]^4)*a[4,1]-6*(a[1,1]^2)*a[3,1]
	dvdc9[27,1]=-6*(a[1,1]^3)*a[4,1]-3*a[1,1]*a[3,1]
	dvdc9[28,1]=-10*(a[1,1]^2)*a[4,1]-a[3,1]
	dvdc9[29,1]=-15*a[1,1]*a[4,1]

	dvdc10[1,1]=0
	dvdc10[2,1]=0
	dvdc10[3,1]=0
	dvdc10[4,1]=0
	dvdc10[5,1]=0
	dvdc10[6,1]=0
	dvdc10[7,1]=0
	dvdc10[8,1]=0
	dvdc10[9,1]=0
	dvdc10[10,1]=0
	dvdc10[11,1]=0
	dvdc10[12,1]=0
	dvdc10[13,1]=0
	dvdc10[14,1]=0
	dvdc10[15,1]=-1
	dvdc10[16,1]=0
	dvdc10[17,1]=-5*a[1,1]*a[2,1]
	dvdc10[18,1]=-(a[2,1]+a[4,1])
	dvdc10[19,1]=0
	dvdc10[20,1]=0
	dvdc10[21,1]=0
	dvdc10[22,1]=0
	dvdc10[23,1]=0
	dvdc10[24,1]=-15*(a[1,1]^2)*a[5,1]
	dvdc10[25,1]=-5*a[1,1]*a[5,1]
	dvdc10[26,1]=-(a[5,1]+a[7,1])
	dvdc10[27,1]=0
	dvdc10[28,1]=0
	dvdc10[29,1]=0

	dvdc11[1,1]=0
	dvdc11[2,1]=0
	dvdc11[3,1]=0
	dvdc11[4,1]=0
	dvdc11[5,1]=0
	dvdc11[6,1]=0
	dvdc11[7,1]=0
	dvdc11[8,1]=0
	dvdc11[9,1]=0
	dvdc11[10,1]=0
	dvdc11[11,1]=0
	dvdc11[12,1]=0
	dvdc11[13,1]=0
	dvdc11[14,1]=0
	dvdc11[15,1]=0
	dvdc11[16,1]=-1
	dvdc11[17,1]=0
	dvdc11[18,1]=0
	dvdc11[19,1]=0
	dvdc11[20,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc11[21,1]=-5*a[1,1]*a[2,1]
	dvdc11[22,1]=0
	dvdc11[23,1]=0
	dvdc11[24,1]=0
	dvdc11[25,1]=0
	dvdc11[26,1]=0
	dvdc11[27,1]=-(a[1,1]^3)*a[5,1]-a[8,1]
	dvdc11[28,1]=-5*(a[1,1]^2)*a[5,1]
	dvdc11[29,1]=-15*a[1,1]*a[5,1]

	dvdc12[1,1]=0
	dvdc12[2,1]=0
	dvdc12[3,1]=0
	dvdc12[4,1]=0
	dvdc12[5,1]=0
	dvdc12[6,1]=0
	dvdc12[7,1]=0
	dvdc12[8,1]=0
	dvdc12[9,1]=0
	dvdc12[10,1]=0
	dvdc12[11,1]=0
	dvdc12[12,1]=0
	dvdc12[13,1]=0
	dvdc12[14,1]=0
	dvdc12[15,1]=0
	dvdc12[16,1]=0
	dvdc12[17,1]=-(a[1,1]^5)
	dvdc12[18,1]=-(a[1,1]^4)
	dvdc12[19,1]=-(a[1,1]^3)
	dvdc12[20,1]=-(a[1,1]^2)
	dvdc12[21,1]=-a[1,1]
	dvdc12[22,1]=0
	dvdc12[23,1]=0
	dvdc12[24,1]=0
	dvdc12[25,1]=0
	dvdc12[26,1]=0
	dvdc12[27,1]=0
	dvdc12[28,1]=0
	dvdc12[29,1]=0

	dvdc13[1,1]=0
	dvdc13[2,1]=0
	dvdc13[3,1]=0
	dvdc13[4,1]=0
	dvdc13[5,1]=0
	dvdc13[6,1]=0
	dvdc13[7,1]=0
	dvdc13[8,1]=0
	dvdc13[9,1]=0
	dvdc13[10,1]=0
	dvdc13[11,1]=0
	dvdc13[12,1]=0
	dvdc13[13,1]=0
	dvdc13[14,1]=0
	dvdc13[15,1]=0
	dvdc13[16,1]=0
	dvdc13[17,1]=0
	dvdc13[18,1]=0
	dvdc13[19,1]=0
	dvdc13[20,1]=0
	dvdc13[21,1]=0
	dvdc13[22,1]=-1
	dvdc13[23,1]=0
	dvdc13[24,1]=-6*a[1,1]*a[2,1]
	dvdc13[25,1]=-(a[2,1]+a[4,1])
	dvdc13[26,1]=0
	dvdc13[27,1]=0
	dvdc13[28,1]=0
	dvdc13[29,1]=0

	dvdc14[1,1]=0
	dvdc14[2,1]=0
	dvdc14[3,1]=0
	dvdc14[4,1]=0
	dvdc14[5,1]=0
	dvdc14[6,1]=0
	dvdc14[7,1]=0
	dvdc14[8,1]=0
	dvdc14[9,1]=0
	dvdc14[10,1]=0
	dvdc14[11,1]=0
	dvdc14[12,1]=0
	dvdc14[13,1]=0
	dvdc14[14,1]=0
	dvdc14[15,1]=0
	dvdc14[16,1]=0
	dvdc14[17,1]=0
	dvdc14[18,1]=0
	dvdc14[19,1]=0
	dvdc14[20,1]=0
	dvdc14[21,1]=0
	dvdc14[22,1]=0
	dvdc14[23,1]=-1
	dvdc14[24,1]=0
	dvdc14[25,1]=0
	dvdc14[26,1]=0
	dvdc14[27,1]=0
	dvdc14[28,1]=-(a[1,1]^2)*a[2,1]-a[3,1]
	dvdc14[29,1]=-6*a[1,1]*a[2,1]

	dvdc15[1,1]=0
	dvdc15[2,1]=0
	dvdc15[3,1]=0
	dvdc15[4,1]=0
	dvdc15[5,1]=0
	dvdc15[6,1]=0
	dvdc15[7,1]=0
	dvdc15[8,1]=0
	dvdc15[9,1]=0
	dvdc15[10,1]=0
	dvdc15[11,1]=0
	dvdc15[12,1]=0
	dvdc15[13,1]=0
	dvdc15[14,1]=0
	dvdc15[15,1]=0
	dvdc15[16,1]=0
	dvdc15[17,1]=0
	dvdc15[18,1]=0
	dvdc15[19,1]=0
	dvdc15[20,1]=0
	dvdc15[21,1]=0
	dvdc15[22,1]=0
	dvdc15[23,1]=0
	dvdc15[24,1]=-(a[1,1]^6)
	dvdc15[25,1]=-(a[1,1]^5)
	dvdc15[26,1]=-(a[1,1]^4)
	dvdc15[27,1]=-(a[1,1]^3)
	dvdc15[28,1]=-(a[1,1]^2)
	dvdc15[29,1]=-a[1,1]

	dvdc=(dvdc1,dvdc2,dvdc3,dvdc4,dvdc5,dvdc6,dvdc7,dvdc8,dvdc9,dvdc10,dvdc11,dvdc12,dvdc13,dvdc14,dvdc15)
	return(dvdc)
}



// Subroutine to compute squeezes.
//
// This subroutine compares the values of the objective function at the
// points c+s*dc and c+0.5*s*dc, with dc = the proposed change in the vector
// of parameters, and with step length s initially at 1. s is halved until
// minus the objective function stops declining.
function squeez4(s_c,s_dc,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
              */ ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
              */ ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
              */ ey3x4_d,ey2x5_d,eyx6_d,maxsqez,obj,estim,neq,w,ID)
{
	s_c1=s_c-s_dc
	s_lm=1/2
	s_itr=1

	if (ID==3 || ID==4)
	{
		s_f1 = EWIDcvec(ey2_d,ex2_d,eyx_d,ey3x_d,ey2x2_d,eyx3_d,ey2x_d,eyx2_d,ID) - s_c1
	}
	else
	{
		s_f1=deff(s_c1,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
			   */ ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
			   */ ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
			   */ ey3x4_d,ey2x5_d,eyx6_d,estim,neq)
	}
	lc1 = s_f1'*w*s_f1

	while (s_itr<=maxsqez)
	{
		s_c2=s_c-s_lm*s_dc
		if (ID==3 || ID==4)
		{
			s_f2 = EWIDcvec(ey2_d,ex2_d,eyx_d,ey3x_d,ey2x2_d,eyx3_d,ey2x_d,eyx2_d,ID) - s_c2
		}
		else
		{
			s_f2=deff(s_c2,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
				   */ ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
				   */ ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
				   */ ey3x4_d,ey2x5_d,eyx6_d,estim,neq)
		}
		lc2 = s_f2'*w*s_f2

		if (lc1 <= lc2 && lc1 <= obj)
		{
			s=(s_c1\s_itr-1)
			return(s)
		}
		else
		{
			s_c1=s_c2
			s_lm=s_lm/2
			lc1=lc2
			s_itr=s_itr+1
		}
	}

	s=(s_c2\s_itr-1)
	return(s)
}

// Subroutine to do the GMM optimuization step
function EWgmm(w,c,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,/*
		    */ ex3_d,ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,/*
			*/ ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,/*
	   	    */ ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,estim,neq,ID)
{
	//-------------set internal optimization variables------------
	maxiter=299
	maxsqez=240
	tol=1e-9

	//-------------start of iteration loop------------------------
	iter=1
	dc=1                     // Initialize the step length.

	while (norm(dc,.)>=tol)
	{
		if (ID==3 || ID==4)
		{
			f = EWIDcvec(ey2_d,ex2_d,eyx_d,ey3x_d,ey2x2_d,eyx3_d,ey2x_d,eyx2_d,ID) - c
			g = -I(neq)
		}
		else
		{
			f=deff(c,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,ey3_d,ey4x_d,/*
				  */ ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,/*
				  */ ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,estim,neq)

			if (estim==2)
			{
				g=grad2(c,neq)
			}
			else if (estim==3)
			{
				g=grad3(c,neq)
			}
			else if (estim==4)
			{
				g=grad4(c,neq)
			}
			else if (estim==5)
			{
				g=grad5(c,neq)
			}
		}

		obj = f'*w*f                // This computes the value of the objective function.

		gwg= g'*w*g                 // This uses the GAUSS-NEWTON method to compute full step dc.
		gwf= g'*w*f
		dc= pinv(gwg)*gwf

		if (maxsqez > 0)
		{
			s=squeez4(c,dc,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
						*/ ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
						*/ ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
						*/ ey3x4_d,ey2x5_d,eyx6_d,maxsqez,obj,estim,neq,w,ID)
			c_new = s[1..rows(c),1]
		}
		else
		{
			c_new=c-dc
		}

		dc=c_new-c                      // Update variables for the next iteration.
		c=c_new

		iter=iter+1
		if (iter>=maxiter)   			// Quit iterating if necessary.
		{
			break
		}
	}

	return(c)
}

function EWIDcvec(ey2_d,ex2_d,eyx_d,ey3x_d,ey2x2_d,eyx3_d,ey2x_d,eyx2_d,ID)
{
	if (ID==3)
	{
		c = J(2,1,0)
		c[1,1] = ey2x_d
		c[2,1] = eyx2_d
	}
	else if (ID==4)
	{
		c = J(8,1,0)
		c[1,1] = ey2_d
		c[2,1] = ex2_d
		c[3,1] = eyx_d
		c[4,1] = ey3x_d
		c[5,1] = ey2x2_d
		c[6,1] = eyx3_d
		c[7,1] = ey2x_d
		c[8,1] = eyx2_d
	}
	else
	{
		errprintf("ID = %d is not a valid option\n", ID)
		exit(error(198))
	}

	return (c)
}

function EWIDoptw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,inezz,zy_d,zx_d,/*
			   */ eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,n,neq,c,ID)
{
	f=J(n,neq,0)

	if (ID==3)
	{
		f[,1] = y2x_d  :- c[1,1] + (-2*eyx_z'*inezz*zy_d'  - ey2_z'*inezz*zx_d')'
		f[,2] = yx2_d  :- c[2,1] + (-ex2_z'*inezz*zy_d'    - 2*eyx_z'*inezz*zx_d')'
	}

	if (ID==4)
	{
		f[,1] = y2_d   :- c[1,1]
		f[,2] = yx_d   :- c[2,1]
		f[,3] = x2_d   :- c[3,1]
		f[,4] = y3x_d  :- c[4,1] + (-3*ey2x_z'*inezz*zy_d' - ey3_z'*inezz*zx_d')'
		f[,5] = y2x2_d :- c[5,1] + (-2*eyx2_z'*inezz*zy_d' - 2*ey2x_z'*inezz*zx_d')'
		f[,6] = yx3_d  :- c[6,1] + (-ex3_z'*inezz*zy_d'    - 3*eyx2_z'*inezz*zx_d')'
		f[,7] = y2x_d  :- c[7,1] + (-2*eyx_z'*inezz*zy_d'  - ey2_z'*inezz*zx_d')'
		f[,8] = yx2_d  :- c[8,1] + (-ex2_z'*inezz*zy_d'    - 2*eyx_z'*inezz*zx_d')'
	}

	return (f)
}


void doEW(string scalar depvar, string scalar misindep, string scalar indepvars, 	///
          string scalar meth, string scalar bXint, 									///
		  string scalar hascons, string scalar touse)
{
	st_view(y, ., depvar, touse)
	st_view(x, ., misindep, touse)
	st_view(z=., ., tokens(indepvars), touse)
	bXint_nums = strtoreal(tokens(bXint))

	n = rows(y)

	if (hascons=="")
	{
		con = J(n,1,1)
		z = (con,z)
	}

	nz = cols(z)

	// Build all the ugly matrices
	{
	mux=pinv(z'*z)*(z'*x)
	muy=pinv(z'*z)*(z'*y)
	y_d=y-z*muy
	x_d=x-z*mux

	y2_d=y_d:^2
	yx_d=(x_d:*y_d)
	x2_d=x_d:^2
	y2x_d=(x_d:*y2_d)
	yx2_d=(x2_d:*y_d)
	y3x_d=y2x_d:*y_d
	y2x2_d=y2_d:*x2_d
	yx3_d=yx2_d:*x_d

	x3_d=x_d:*x2_d
	y3_d=y_d:*y2_d
	y4x_d=y3_d:*yx_d
	y3x2_d=y3_d:*x2_d
	y2x3_d=y2_d:*x3_d
	yx4_d=yx_d:*x3_d
	y4_d=y2_d:^2
	x4_d=x2_d:^2
	y5x_d=y_d:*y4x_d
	y4x2_d=y_d:*y3x2_d
	y3x3_d=y_d:*y2x3_d
	y2x4_d=y_d:*yx4_d
	yx5_d=yx4_d:*x_d

	y5_d=y2_d:*y3_d
	x5_d=x2_d:*x3_d
	y6x_d=y5_d:*yx_d
	y5x2_d=y5_d:*x2_d
	y4x3_d=y4_d:*x3_d
	y3x4_d=y3_d:*x4_d
	y2x5_d=y2_d:*x5_d
	yx6_d=yx_d:*x5_d

	ey2_d=mean(y2_d)
	eyx_d=mean(yx_d)
	ex2_d=mean(x2_d)
	ey2x_d=mean(y2x_d)
	eyx2_d=mean(yx2_d)
	ey3x_d=mean(y3x_d)
	ey2x2_d=mean(y2x2_d)
	eyx3_d=mean(yx3_d)

	ex3_d=mean(x3_d)
	ey3_d=mean(y3_d)
	ey4x_d=mean(y4x_d)
	ey3x2_d=mean(y3x2_d)
	ey2x3_d=mean(y2x3_d)
	eyx4_d=mean(yx4_d)
	ey4_d=mean(y4_d)
	ex4_d=mean(x4_d)
	ey5x_d=mean(y5x_d)
	ey4x2_d=mean(y4x2_d)
	ey3x3_d=mean(y3x3_d)
	ey2x4_d=mean(y2x4_d)
	eyx5_d=mean(yx5_d)

	ey5_d=mean(y5_d)
	ex5_d=mean(x5_d)
	ey6x_d=mean(y6x_d)
	ey5x2_d=mean(y5x2_d)
	ey4x3_d=mean(y4x3_d)
	ey3x4_d=mean(y3x4_d)
	ey2x5_d=mean(y2x5_d)
	eyx6_d=mean(yx6_d)


	// Make the moments for the correct standard errors for beta.

	y6_d=y_d:^6
	x6_d=x_d:^6

	idx = J(nz,1,1)
	ey2_z   = mean(y2_d[,idx]:*z)'
	eyx_z   = mean(yx_d[,idx]:*z)'
	ex2_z   = mean(x2_d[,idx]:*z)'
	ey2x_z  = mean(y2x_d[,idx]:*z)'
	eyx2_z  = mean(yx2_d[,idx]:*z)'
	ey3x_z  = mean(y3x_d[,idx]:*z)'
	ey2x2_z = mean(y2x2_d[,idx]:*z)'
	eyx3_z  = mean(yx3_d[,idx]:*z)'
	ey3_z   = mean(y3_d[,idx]:*z)'
	ex3_z   = mean(x3_d[,idx]:*z)'
	ey4x_z  = mean(y4x_d[,idx]:*z)'
	ey3x2_z = mean(y3x2_d[,idx]:*z)'
	ey2x3_z = mean(y2x3_d[,idx]:*z)'
	eyx4_z  = mean(yx4_d[,idx]:*z)'
	ey4_z   = mean(y4_d[,idx]:*z)'
	ex4_z   = mean(x4_d[,idx]:*z)'
	ey5x_z  = mean(y5x_d[,idx]:*z)'
	ey4x2_z = mean(y4x2_d[,idx]:*z)'
	ey3x3_z = mean(y3x3_d[,idx]:*z)'
	ey2x4_z = mean(y2x4_d[,idx]:*z)'
	eyx5_z  = mean(yx5_d[,idx]:*z)'
	ey5_z   = mean(y5_d[,idx]:*z)'
	ex5_z   = mean(x5_d[,idx]:*z)'
	ey6x_z  = mean(y6x_d[,idx]:*z)'
	ey5x2_z = mean(y5x2_d[,idx]:*z)'
	ey4x3_z = mean(y4x3_d[,idx]:*z)'
	ey3x4_z = mean(y3x4_d[,idx]:*z)'
	ey2x5_z = mean(y2x5_d[,idx]:*z)'
	eyx6_z  = mean(yx6_d[,idx]:*z)'
	ey6_z   = mean(y6_d[,idx]:*z)'
	ex6_z   = mean(x6_d[,idx]:*z)'

	ezz = (z'*z)/n
	inezz = pinv(ezz)
	zy_d = z:*y_d[,idx]
	zx_d = z:*x_d[,idx]
	}

	//----------here we do identification tests.------------------
	for (ID=3;ID<=4;ID++)
	{
		c   = EWIDcvec(ey2_d,ex2_d,eyx_d,ey3x_d,ey2x2_d,eyx3_d,ey2x_d,eyx2_d,ID)
		neq = rows(c)
		ff  = EWIDoptw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,inezz,zy_d,zx_d,/*
			        */ eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,n,neq,c,ID)
		w   = pinv((ff'*ff)/n)
		c   = EWgmm(w,c,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,/*
		         */ ex3_d,ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,/*
				 */ ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,/*
				 */ ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,-1,neq,ID)
		ff  = EWIDoptw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,inezz,zy_d,zx_d,/*
			        */ eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,n,neq,c,ID)
		w   = pinv((ff'*ff)/n)
		if (ID==3)
		{
			ID3val  = n*c'*w*c     		//Newey West Wald test for c==0
			ID3stat = 1-chi2(2,ID3val)
		}
		if (ID==4)
		{
			Cum4 = J(5,1,0)
			Cum4[1,1] = ey3x_d  - 3*eyx_d*ey2_d
			Cum4[2,1] = ey2x2_d - ey2_d*ex2_d - 2*eyx_d^2
			Cum4[3,1] = eyx3_d  - 3*eyx_d*ex2_d
			Cum4[4,1] = ey2x_d
			Cum4[5,1] = eyx2_d

			Vee = J(5,neq,0)
			Vee[1,1] = -3*eyx_d
			Vee[1,2] = -3*ey2_d
			Vee[1,4] = 1
			Vee[2,1] = -ex2_d
			Vee[2,2] = -4*eyx_d
			Vee[2,3] = -ey2_d
			Vee[2,5] = 1
			Vee[3,2] = -3*ex2_d
			Vee[3,3] = -3*eyx_d
			Vee[3,6] = 1
			Vee[4,7] = 1
			Vee[5,8] = 1

			w = pinv(Vee*((ff'*ff)/n)*Vee')
			ID4val  = n*Cum4'*w*Cum4
			ID4stat = 1-chi2(5,ID4val)
		}
	}


	// Now define method

	if (strupper(meth)=="GMM3" || meth=="")
	{
		nc=5
		neq=5
		estim=1
	}
	else if (strupper(meth)=="GMM4")
	{
		nc=6
		neq=8
		estim=2
	}
	else if (strupper(meth)=="GMM5")
	{
		nc=9
		neq=14
		estim=3
	}
	else if (strupper(meth)=="GMM6")
	{
		nc=12
		neq=21
		estim=4
	}
	else if (strupper(meth)=="GMM7")
	{
		nc=15
		neq=29
		estim=5
	}
	else
	{
		errprintf("method = %s is not a valid option\n", meth)
		exit(error(198))
	}

	// Start the GMM part
	//----------here we input the starting values.----------------

	if (estim==1)
	{
		bXint_nums = 0
	}

	c       = J(nc,1,0)
	bXlen   = cols(bXint_nums)
	objsave = J(bXlen,1,0)

	for (rep=1;rep<=(bXlen+1);rep++)
	{
		if (rep<(bXlen+1))
		{
			c[1,1] = bXint_nums[1,rep]
		}
		else
		{
			if (bXlen>1)
			{
				// restore best inital value
				minindex(objsave, 1, ind, where)
				c[1,1] = bXint_nums[1,ind[1,1]]
			}
			else
			{
				// no point in restoring - only one value tested
				break
			}
		}

		if (c[1,1]==0)
		{
			c[1,1]=(ey2x_d)/(eyx2_d)
		}

		c[2,1]=(eyx_d)/c[1,1]
		c[3,1]=(ey2_d)-(c[1,1]^2)*c[2,1]
		c[4,1]=(ex2_d)-c[2,1]
		c[5,1]=(eyx2_d)/c[1,1]

		if (estim>1)
		{
			c[6,1]=((eyx3_d)-3*c[1,1]*c[2,1]*c[4,1])/c[1,1]
		}
		if (estim>2)
		{
			c[7,1]=(ex3_d)-c[5,1]
			c[8,1]=(ey3_d)-(c[1,1]^3)*c[5,1]
			c[9,1]=((eyx4_d)/c[1,1])-6*c[5,1]*c[4,1]-4*c[2,1]*c[7,1]
		}
		if (estim>3)
		{
			c[10,1]=(ey4_d)-(c[1,1]^4)*c[6,1]-6*(c[1,1]^2)*c[2,1]*c[3,1]
			c[11,1]=(ex4_d)-c[6,1]-6*c[2,1]*c[4,1]
			c[12,1]=((eyx5_d)/c[1,1])-10*c[6,1]*c[4,1]-10*c[5,1]*c[7,1] ///
					-5*c[2,1]*c[11,1]
		}
		if (estim>4)
		{
			c[13,1]=(ey5_d)-(c[1,1]^5)*c[9,1]-10*(c[1,1]^3)*c[5,1]*c[3,1] ///
					-10*(c[1,1]^2)*c[2,1]*c[8,1]
			c[14,1]=(ex5_d)-c[9,1]-10*c[5,1]*c[4,1]-10*c[2,1]*c[7,1]
			c[15,1]=((eyx6_d)/c[1,1])-15*c[9,1]*c[4,1]-20*c[6,1]*c[7,1] ///
					-15*c[5,1]*c[11,1]-6*c[2,1]*c[14,1]
		}

		// We only need to maximize if the estimator is not exactly identified -> not
		// for GMM3, but for all the rest. Note that gradX calculates analytical
		// derivetives
		if (estim>1)
		{
			//--------the program creates the weighting matrix.------------
			ff=optw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,x3_d,y3_d,y4x_d,/*
				  */y3x2_d,y2x3_d,yx4_d,y4_d,x4_d,y5x_d,y4x2_d,y3x3_d,y2x4_d,yx5_d,/*
				  */y5_d,x5_d,y6x_d,y5x2_d,y4x3_d,y3x4_d,y2x5_d,yx6_d,/*
				  */ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
				  */ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
				  */ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
				  */ey3x4_d,ey2x5_d,eyx6_d,/*
				  */inezz,zy_d,zx_d,eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,ey4_z,/*
				  */ey3x_z,ey2x2_z,eyx3_z,ex4_z,ey5_z,ey4x_z,ey3x2_z,ey2x3_z,eyx4_z,ex5_z,/*
				  */ey6_z,ey5x_z,ey4x2_z,ey3x3_z,ey2x4_z,eyx5_z,ex6_z,n,neq,estim)
			w=pinv((ff'*ff)/n)

			c = EWgmm(w,c,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,/*
				   */ ex3_d,ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,/*
				   */ ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,/*
				   */ ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,estim,neq,0)
		}

		// Now do t-stats and friends

		ff=optw(y2_d,yx_d,x2_d,y2x_d,yx2_d,y3x_d,y2x2_d,yx3_d,x3_d,y3_d,y4x_d,/*
			  */y3x2_d,y2x3_d,yx4_d,y4_d,x4_d,y5x_d,y4x2_d,y3x3_d,y2x4_d,yx5_d,/*
			  */y5_d,x5_d,y6x_d,y5x2_d,y4x3_d,y3x4_d,y2x5_d,yx6_d,/*
			  */ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,/*
			  */ey3_d,ey4x_d,ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,/*
			  */ey3x3_d,ey2x4_d,eyx5_d,ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,/*
			  */ey3x4_d,ey2x5_d,eyx6_d,/*
			  */inezz,zy_d,zx_d,eyx_z,ey2_z,ex2_z,ey3_z,ey2x_z,eyx2_z,ex3_z,ey4_z,/*
			  */ey3x_z,ey2x2_z,eyx3_z,ex4_z,ey5_z,ey4x_z,ey3x2_z,ey2x3_z,eyx4_z,ex5_z,/*
			  */ey6_z,ey5x_z,ey4x2_z,ey3x3_z,ey2x4_z,eyx5_z,ex6_z,n,neq,estim)
		f=deff(c,ey2_d,eyx_d,ex2_d,ey2x_d,eyx2_d,ey3x_d,ey2x2_d,eyx3_d,ex3_d,ey3_d,ey4x_d,/*
			  */ ey3x2_d,ey2x3_d,eyx4_d,ey4_d,ex4_d,ey5x_d,ey4x2_d,ey3x3_d,ey2x4_d,eyx5_d,/*
			  */ ey5_d,ex5_d,ey6x_d,ey5x2_d,ey4x3_d,ey3x4_d,ey2x5_d,eyx6_d,estim,neq)
		w=pinv((ff'*ff)/n)
		obj = f'*w*f

		if (rep<(bXlen+1))
		{
			objsave[rep,1]=obj
		}
	}


	if (estim==1)
	{
		g=grad1(c,neq)
	}
	else if (estim==2)
	{
		g=grad2(c,neq)
	}
	else if (estim==3)
	{
		g=grad3(c,neq)
	}
	else if (estim==4)
	{
		g=grad4(c,neq)
	}
	else if (estim==5)
	{
		g=grad5(c,neq)
	}

	gwg=g'*w*g
	vcX=pinv(gwg)/n
	stderr=sqrt(diagonal(vcX))
	inflnc=-n*vcX*g'*w*ff'

	bX  = c[1,1]
	seX = stderr[1,1]
	inflncX = inflnc[1,]'
	dfree = neq-nc

	if (estim>1)
	{
		Jstat = obj*n
		Jval  = 1-chi2(dfree,Jstat)
	}
	else
	{
		Jstat = -1
		Jval  = -1
	}

	// Make bZ, seZ and friends

	c11 = c[1,1]
	bZ = muy[,1]-c11*mux[,1]

	gee = J((nz*2+1), nz, 0)
	gee[1..nz,1..nz] = I(nz)
	gee[(nz+1)..(nz*2),1..nz] = -c11*I(nz)
	gee[(nz*2+1),] = -mux[,1]'

	bigphi = ((inezz*zy_d')\(inezz*zx_d')\(-inflnc[1,]))
	avar = (bigphi*bigphi')/n:^2
	vcZ = gee'*avar*gee
	seZ = sqrt(diagonal(vcZ))
	inflncZ = -bigphi'*gee

	//--------------rho2 and tau2-------------------------------------------------
	ez = mean(z)
	sigz=(z-J(n,1,1)#ez)'*(z-J(n,1,1)#ez)/n

	numer=muy'*sigz*muy + c[1,1]:^2*c[2,1]
	rho=numer/(numer+c[3,1])

	numer=mux'*sigz*mux + c[2,1]
	tau=numer/(numer+c[4,1])

	vecsigz=J(nz+nz*(nz-1)/2,1,0)
	vecsigz[1..nz,1]=diagonal(sigz)

	phiz=J(n,nz+nz*(nz-1)/2,0)
	phiz[,1..nz]=(z-J(n,1,1)#ez):*(z-J(n,1,1)#ez)

	// here we pack the top half of the sigz matrix into a vector
	counter=nz+nz*(nz-1)/2
	for (qq=nz;qq>=2;qq--)
	{
		for (ee=qq-1;ee>=1;ee--)
		{
			vecsigz[counter,1]=sigz[ee,qq]
			phiz[,counter]=(z[,ee]-J(n,1,1)#ez[1,ee]):*(z[,qq]-J(n,1,1)#ez[1,qq])
			counter=counter-1
		}
	}

	phiz=phiz-J(n,1,1)#vecsigz'
	phimuy=inezz*zy_d'
	phimux=inezz*zx_d'

	//--make the influence functions for the standard errors for rho2 and tau2.
	bigphi=(phimux\phimuy\(phiz')\(-inflnc[1..4,]))
	avar=(bigphi*bigphi')/n:^2
	//--------------first column is for rho2 and the second is for tau2.
	gee = J(2*nz+(nz+nz*(nz-1)/2)+4,2,0)
	//-----------first,do rho2---------------------------
	numer=muy'*sigz*muy + c[1,1]:^2*c[2,1]
	denom=numer + c[3,1]; ;
	//------------derivatives wrt muy---------------------
	for (qq=1;qq<=nz;qq++)
	{
		gee[nz+qq,1]= (2*muy'*sigz[,qq])/denom - numer*(2*muy'*sigz[,qq])/(denom:^2)
	}
	//------------derivatives wrt the first part of sigz--
	for (qq=1;qq<=nz;qq++)
	{
		gee[2*nz+qq,1]=(muy[qq,1]:^2)/denom-numer/(denom:^2)*muy[qq,1]:^2
	}
	//------------derivatives wrt the second part of sigz--
	counter=2*nz+(nz+nz*(nz-1)/2)
	for (qq=nz;qq>=2;qq--)
	{
		for (ee=qq-1;ee>=1;ee--)
		{
			gee[counter,1]= 2*muy[ee,1]*muy[qq,1]/denom-2*numer/(denom:^2)*muy[ee,1]*muy[qq,1]
			counter=counter-1
		}
	}
	//------------derivatives wrt c--
	counter=2*nz+(nz+nz*(nz-1)/2)
	gee[counter+1,1]=2*c[1,1]*c[2,1]/denom-2*numer/(denom:^2)*c[1,1]*c[2,1]
	gee[counter+2,1]=(c[1,1]:^2)/denom-numer/(denom:^2)*c[1,1]:^2
	gee[counter+3,1]=-numer/(denom:^2)

	//-----------now for tau2------------------------
	numer=mux'*sigz*mux+c[2,1]
	denom=numer+c[4,1]
	//------------derivatives wrt mux---------------------
	for (qq=1;qq<=nz;qq++)
	{
		gee[qq,2]=(2*mux'*sigz[,qq])/denom-numer*(2*mux'*sigz[,qq])/(denom:^2)
	}
	//------------derivatives wrt the first part of sigz--
	for (qq=1;qq<=nz;qq++)
	{
		gee[2*nz+qq,2]=(mux[qq,1]:^2)/denom-numer/(denom:^2)*mux[qq,1]:^2
	}
	//------------derivatives wrt the second part of sigz--
	counter=2*nz+(nz+nz*(nz-1)/2)
	for (qq=nz;qq>=2;qq--)
	{
		for (ee=qq-1;ee>=1;ee--)
		{
			gee[counter,2]= 2*mux[ee,1]*mux[qq,1]/denom-2*numer/(denom:^2)*mux[ee,1]*mux[qq,1]
			counter=counter-1
		}
	}
	//------------derivatives wrt c--
	counter=2*nz+(nz+nz*(nz-1)/2)
	gee[counter+2,2]=1/denom-numer/(denom:^2)
	gee[counter+4,2]=-numer/(denom:^2)


	//--------------done with influnence funcs. Now use them.----------------------

	vcrhotau=gee'*avar*gee

	inflncrho=-bigphi'*gee[,1]
	inflnctau=-bigphi'*gee[,2]

	SErho=sqrt(vcrhotau[1,1])
	SEtau=sqrt(vcrhotau[2,2])

	//--------------All done. Save return values and exit.--------------------------
	beta = J(nz+1,1,1)
	beta[1,1] = bX
	beta[2..(nz+1),1] = bZ
	inflnc = inflncX, inflncZ
	vcmat = inflnc'*inflnc/n:^2
	serr = sqrt(diagonal(vcmat))
	st_matrix("r(beta)",beta')
	st_matrix("r(VCmat)",vcmat)
	st_matrix("serr",serr)
	st_matrix("inflnc",inflnc)
	st_numscalar("r(N)",n)
	st_numscalar("r(Jstat)",Jstat)
	st_numscalar("r(Jval)",Jval)
	st_numscalar("r(dfree)",dfree)
	st_numscalar("r(rho)",rho)
	st_numscalar("r(tau)",tau)
	st_numscalar("r(SErho)",SErho)
	st_numscalar("r(SEtau)",SEtau)
	st_matrix("vcrhotau",vcrhotau)
	st_matrix("inflncrhotau",(inflncrho, inflnctau))
	st_matrix("w",w)
	st_numscalar("r(obj)",obj)
	st_numscalar("r(ID3val)",ID3val)
	st_numscalar("r(ID3stat)",ID3stat)
	st_numscalar("r(ID4val)",ID4val)
	st_numscalar("r(ID4stat)",ID4stat)
}
end