/*******************************************************************************
** (C) KEISUKE KONDO
**
** Release Date: November 15, 2016
** Update Date: July 11, 2022
** Version: 1.20
**
** [Reference]
** Kondo, K. (2017) "Quantile approach for distinguishing agglomeration from
** firm selection in Stata," RIETI TP 17-T-001.
** GitHub: https://github.com/keisukekondokk/estquant
**
** Combes, P.P., Duranton, G., Gobillon, L., Puga, D., and Roux, S., (2012)
** "The Productivity Advantages of Large Cities: Distinguishing Agglomeration
** From Firm Selection," Econometrica 80(6), pp. 2543-2594
** URL: https://doi.org/10.3982/ECTA8442
**
** [Contact]
** Email: kondo-keisuke@rieti.go.jp
** URL: https://keisukekondokk.github.io/
********************************************************************************/
** Version: 1.20
** Improved program code
** Added OPTECHnique() option
** Added OPNMSimplexdeltas() option
** Version: 1.10
** Added INITR() option
** Version: 1.04
** Improved program code
** Version: 1.03
** Bug fix (Delete "mata mlib create")
** Version: 1.02
** Integrated program code into one ado file
** Version: 1.01
** Improved program code
**
**

/*
** ==========================================================
** MATA
** ==========================================================
*/
clear mata
/*---------- START Class quantilef ----------*/
version 12
mata:
mata set matastrict off
class quantilef {
	/* Member Variables */
	real scalar truncation, shift, dilation
	real scalar cref1, cref2, nrange
	string scalar optechnique
	real scalar opnmsdeltas
	real scalar maxIt, eps1, eps2
	real matrix mP
	real vector vbeta, beta0, betaf
	real scalar pgammf
	real vector y1, y2
	real scalar mean_y1, mean_y2, mean_y
	real scalar sd_y1, sd_y2, sd_y
	real matrix y
	real scalar Cr, Ct1, Ct2
	real scalar nobs, nobs_c1, nobs_c2
	real scalar myn1, myn2, sdyn

	/* Member Functions */
	void new()
	void init_setting_range(), init_setting_data(), init_setting_category(), /*
		*/ init_setting_spec(), init_setting_params(), init_setting_conv()
	real matrix zb(), dzb_b(), zb_1(), dzb_1_b()
	real matrix zb0(), zb0_1()
	real scalar dzb_u()
	real scalar gammapp(), gammapp_1()
	real matrix tg(), dtg_g(), tg_1(), dtg_u()
	real matrix weight1(), weight2()
	real matrix quant()
	real scalar init()
	void critb()
	void critbisf()
	real scalar critb0(), critbisf0()
	real matrix BlGart()
	real scalar calctest()
}
/* =====================================
  Member Functions
===================================== */
/* Constructor */
void quantilef::new()
{
}
/* Initial Setting Range */
void quantilef::init_setting_range( real scalar vRange )
{
	nrange = vRange
}
/* Initial Setting Data */
void quantilef::init_setting_data( real vector vY1, real vector vY2 )
{
	y1 = vY1
	y2 = vY2
	mean_y1 = mean(y1)
	mean_y2 = mean(y2)
	mean_y = mean(y1 \ y2)
	sd_y1 = sqrt(variance(y1))
	sd_y2 = sqrt(variance(y2))
	sd_y = sqrt(variance(y1 \ y2))
	nobs_c1 = rows(y1)
	nobs_c2 = rows(y2)
	nobs = nobs_c1 + nobs_c2
}
/* Initial Setting Category */
void quantilef::init_setting_category( numeric vector vCat )
{
	cref1 = vCat[1]
	cref2 = vCat[2]
}
/* Initial Setting Specification */
void quantilef::init_setting_spec( numeric vector vSpec )
{
	truncation = vSpec[1]
	shift = vSpec[2]
	dilation = vSpec[3]
}
/* Initial Setting Parameter */
void quantilef::init_setting_params( real rowvector vBeta )
{
	/* vbeta = (D, A) */
	beta0 = vBeta
	betaf = vBeta
	vbeta = J(1, 2, 1)
	vbeta[1] = 1
	vbeta[2] = 0
}
/* Initial Setting for Convergnece */
void quantilef::init_setting_conv( real scalar maxIt, real scalar eps1, real scalar eps2 )
{
	maxIt = maxIt
	eps1 = eps1
	eps2 = eps2
}
/*
** ==========================================================
** 1. Specify Parametric Functions of Value
** ==========================================================
*/
/* Value Transformation Function */
real matrix quantilef::zb( real vector u, real vector beta )
{
	real matrix zb
	if( dilation == 1 & shift == 0 ){
		zb = beta :* u
	}
	else if( dilation == 0 & shift == 1 ){
		zb = u :+ beta
	}
	else if( dilation == 1 & shift == 1 ){
		zb = beta[1] :* u :+ beta[2]
	}
	else if( dilation == 0 & shift == 0 ){
		zb = u
	}
	return(zb)
}
real matrix quantilef::zb0( real vector u )
{
	/* (D, A) = (1, 0) */
	real matrix zb0
	zb0 = u
	return(zb0)
}
/* Partial Derivative of zb() wrt beta */
real matrix quantilef::dzb_b( real vector u, real vector beta )
{
	real matrix dzb
	if( dilation == 1 & shift == 0 ){
		dzb = u
	}
	else if( dilation == 0 & shift == 1 ){
		dzb = J(rows(u), 1, 1)
	}
	else if( dilation == 1 & shift == 1 ){
		dzb = ( u, J(rows(u), 1, 1) )
	}
	return(dzb)
}
/* Partial Derivative of zb() wrt u */
real scalar quantilef::dzb_u( real vector u, real vector beta )
{
	real scalar dzb
	if( dilation == 1 & shift == 0 ){
		dzb = beta
	}
	else if( dilation == 0 & shift == 1 ){
		dzb = 1
	}
	else if( dilation == 1 & shift == 1 ) {
		dzb = beta[1]
	}
	return(dzb)
}
/* Inverse of zb(), Holding beta Constant */
real matrix quantilef::zb_1( real vector x, real vector beta )
{
	real matrix zb_1
	if( dilation == 1 & shift == 0 ){
		zb_1 = x :/ beta
	}
	else if( dilation == 0 & shift == 1 ){
		zb_1 = x :- beta
	}
	else if( dilation == 1 & shift == 1 ) {
		zb_1 = ( x :- beta[2] ) :/ beta[1]
	}
	else if( dilation == 0 & shift == 0 ) {
		zb_1 = x
	}
	return(zb_1)
}
real matrix quantilef::zb0_1( real vector x )
{
	/* (D, A) = (1, 0) */
	real matrix zb0_1
	zb0_1 = x
	return(zb0_1)
}
/* Derivative of zb_1() wrt beta */
real matrix quantilef::dzb_1_b( real vector x, real vector beta )
{
	real matrix dzb
	if( dilation == 1 & shift == 0 ){
		dzb = - x :/ (beta^2)
	}
	else if( dilation == 0 & shift == 1 ){
		dzb = - J(rows(x), 1, 1)
	}
	else if( dilation == 1 & shift == 1 ){
		dzb = ( -(x :- beta[2]) :/ (beta[1]^2), -(1 :/ beta[1]) :* J(rows(x), 1, 1) )
	}
	return(dzb)
}
/*
** ==========================================================
** 2. Specify Functions of Rank Transformation or Selection
** =========================================================
*/
/* Selection */
real scalar quantilef::gammapp( real scalar pgamm )
{
	real scalar gamm
	gamm = ( exp(pgamm) - 5 ) / ( 5 + exp(pgamm) )
	gam = gamm
	return(gam)
}
/* Selection */
real scalar quantilef::gammapp_1( real scalar gamm )
{
	real scalar pgamm
	gam = gamm
	if( gam == 1 ){
		/* S = 1 */
		gam = gamm - (1e-8)
	}
	pgamm = log( (1 + gam) / (1 - gam) ) + log(5)
	return(pgamm)
}
/* Selection */
real matrix quantilef::tg( real vector u, real scalar gamma )
{
	real matrix tg
	gam = gamma
	tg = gam :+ (1 - gam) :* u
	return(tg)
}
/* Selection */
real matrix quantilef::dtg_g( real vector u, real scalar gamma)
{
	real matrix dtg
	dtg = 1 :- u
	return(dtg)
}
/* Selection */
real matrix quantilef::dtg_u(real vector u, real scalar gamma)
{
	real matrix dtg
	gam = gamma
	dtg = J(rows(u), 1, (1 - gam))
	return(dtg)
}
/* Selection */
real matrix quantilef::tg_1(real vector x, real scalar gamma)
{
	real matrix tg_1
	gam = gamma
	if( gam == 1 ){
		/* S = 1 */
		gam = gamma - (1e-8)
	}
	tg_1 = ( x :- gam ) :/ ( 1 - gam )
	return(tg_1)
}
/*
** ==========================================================
** Weight and Quantile
** ==========================================================
*/
/* Weight 1 */
real matrix quantilef::weight1(real vector x)
{
	real matrix w
	w = J(rows(x), cols(x), 1)
	return(w)
}
/* Weight 2 */
real matrix quantilef::weight2(real vector x)
{
	real matrix w
	w = J(rows(x), cols(x), 1)
	return(w)
}
/* Procedure estimating quantiles of rank (rang2) from the distribution of values in y1 using linear interpolation */
real matrix quantilef::quant(real vector vy1, real vector rang2)
{
	real matrix qy1
	nobs1 = rows(vy1)
	nobs2 = rows(rang2)
	intb = trunc( nobs1 :* rang2 )
	intb = intb :+ (intb :== J(nobs2, 1, 0))
	intu = trunc( nobs1 :* rang2 ) :+ 1
	intu = intu :- (intu :> J(nobs2, 1, nobs1))
	qy1b = vy1[intb]
	qy1u = vy1[intu]
	qy1 = ( ((nobs1 :* rang2) :- intb) :* (qy1u :- qy1b) ) :+ qy1b
	return(qy1)
}
/*
** ==========================================================
** Optimize Criteria wrt Beta (not Gamma)
** ==========================================================
*/
/* Optimizing the Criteria wrt beta (and not gamma) */
real scalar quantilef::init(class quantilef Q)
{
	if( dilation == 1 & shift == 0 ){
		vbeta[1] = sdyn
		vbeta[2] = 0
		beta0 = vbeta[1]
	}
	else if( dilation == 0 & shift == 1 ) {
		vbeta[1] = 1
		vbeta[2] = myn2 - myn1
		beta0 = vbeta[2]
	}
	else {
		vbeta[1] = sdyn
		vbeta[2] = myn2 - myn1 * sdyn
		beta0 = vbeta
	}

	/* Optimization */
	Sb = optimize_init()
	optimize_init_which(Sb, "min")
	optimize_init_evaluator(Sb, &eval_critb())
	optimize_init_evaluatortype(Sb, "d1")
	optimize_init_argument(Sb, 1, Q)
	if(optechnique == "nr" | optechnique == "dfp" | optechnique == "bfgs"){
		optimize_init_technique(Sb, optechnique)
		optimize_init_singularHmethod(Sb, "hybrid")
	}
	else if(optechnique == "nm") {
		optimize_init_technique(Sb, optechnique)
		optimize_init_nmsimplexdeltas(Sb, opnmsdeltas)
	}
	optimize_init_params(Sb, beta0)
	optimize_init_tracelevel(Sb, "none")
	optimize_init_conv_maxiter(Sb, maxIt)
	optimize_init_conv_ptol(Sb, 1e-20)
	optimize_init_conv_vtol(Sb, eps1)
	optimize_init_conv_nrtol(Sb, eps2)
	betaf = optimize(Sb)
	cr = optimize_result_value(Sb)
	return(cr)
}
/*
** ==========================================================
** Estimating Parameters for the considered transformation
** ==========================================================
*/
/* Computation of the Criterium to Minimize, as a function of beta only */
void quantilef::critb(todo, beta, v, g, H)
{
	posb = range(1, (rows(y) - 1), 1)
	posu = range(2, rows(y), 1)
	ys = sort(y, 3)
	w1 = ys[., 5]
	yt1 = ys[., 1]
	yt2 = ys[., 2]
	rgt1 = ys[., 3]
	It1 = (w1 :* (zb(yt1, beta) :- yt2)):^2
	Ct1 = ( (It1[posu] :+ It1[posb])'*(rgt1[posu] :- rgt1[posb]) ) / (rgt1[rows(y)] - rgt1[1]) / 2
	dIt1 = 2 :* ( J(1, length(beta), 1)#(w1 :* (zb(yt1, beta) :- yt2)) ) :* dzb_b(yt1, beta)
	dCt1 = ( (dIt1[posu, .] :+ dIt1[posb, .])'*(rgt1[posu] :- rgt1[posb]) ) / (rgt1[rows(y)] - rgt1[1]) / 2
	ys = sort(y, 4)
	w2 = ys[., 6]
	yt1 = ys[., 1]
	yt2 = ys[., 2]
	rgt2 = ys[., 4]
	It2 = (w2 :* (yt1 :- zb_1(yt2, beta))):^2
	Ct2 = ( (It2[posu] :+ It2[posb])'*(rgt2[posu] :- rgt2[posb]) ) / (rgt2[rows(y)] - rgt2[1]) / 2
	dIt2 = -2 :* ( J(1, length(beta), 1)#(w2 :* (yt1 :- zb_1(yt2, beta))) ) :* dzb_1_b(yt2, beta)
	dCt2 = ( (dIt2[posu, .] :+ dIt2[posb, .])'*(rgt2[posu] :- rgt2[posb]) ) / (rgt2[rows(y)] - rgt2[1]) / 2
	v = Ct1 + Ct2
	/* Derivative of the Criterium */
	if( todo >= 1 ){
		g = (dCt1 :+ dCt2)'
	}
}
/* Computation of the Criterium as a Function of pgamm */
void quantilef::critbisf(todo, p, Q, v, g, H)
{
	/* Get Gamma */
	pgamm = p
	gamm = gammapp(pgamm)
	if( truncation == 0 ){
		gamm = 0
	}

	/* common support condition */
	umin1 = max( (0, tg_1(0, gamm)) )
	umax1 = min( (1, tg_1(1, gamm)) )
	umin2 = max( (0, tg(0, gamm)) )
	umax2 = min( (1, tg(1, gamm)) )

	/* rank transformation */
	n1t = trunc( rows(y1) :* (umax1 - umin1) )
	n2t = trunc( rows(y2) :* (umax2 - umin2) )

	/* computation of threshold ranks */
	rangn1 = umin1 :+ ((range(0, nrange, 1)) :/ nrange) :* (umax1 :- umin1)
	tgrangn1 = tg( rangn1, gamm )
	rangn2 = umin2 :+ ((range(0, nrange, 1)) :/ nrange) :* (umax2 :- umin2)
	tg_rangn2 = tg_1( rangn2, gamm )
	w1 = weight1( rangn1 \ tg_rangn2 )
	w2 = weight2( tgrangn1 \ rangn2 )
	yn1 = quant( y1, rangn1 )
	qy2 = quant( y2, tgrangn1 )
	yn2 = quant( y2, rangn2 )
	qy1 = quant( y1, tg_rangn2 )
	y = ((yn1, qy2, rangn1, tgrangn1) \ (qy1, yn2, tg_rangn2, rangn2)), /*
		*/ ( w1, w2, (J(rows(yn1), 1, cref1) \ J(rows(yn2), 1, cref2)) )

	/* Optimization wrt beta (given gamma) */
	if( dilation == 0 & shift == 0 ){
		Cr = critb0()
	}
	else if( dilation == 1 | shift == 1 ){
		/* Calculate Initial Values for Optimization */
		myn1 = mean(yn1)
		myn2 = mean(yn2)
		vyn1 = variance(yn1)
		vyn2 = variance(yn2)
		sdyn = 1
		if( vyn1 > 0 ){
			sdyn = sqrt( vyn2 / vyn1 )
		}
		/* Optimization wrt beta (given gamma) */
		Cr = init(Q)
	}

	/* Criterium */
	v = Cr
}
/* Computation of the Criterium for R2 */
real scalar quantilef::critb0()
{
	posb = range(1, (rows(y) - 1), 1)
	posu = range(2, rows(y), 1)
	ys = sort(y, 3)
	w1 = ys[., 5]
	yt1 = ys[., 1]
	yt2 = ys[., 2]
	rgt1 = ys[., 3]
	It1 = (w1 :* (zb0(yt1) :- yt2)):^2
	Ct1 = ( (It1[posu] :+ It1[posb])'*(rgt1[posu] :- rgt1[posb]) ) / (rgt1[rows(y)] - rgt1[1]) / 2
	ys = sort(y, 4)
	w2 = ys[., 6]
	yt1 = ys[., 1]
	yt2 = ys[., 2]
	rgt2 = ys[., 4]
	It2 = (w2 :* (yt1 :- zb0_1(yt2))):^2
	Ct2 = ( (It2[posu] :+ It2[posb])'*(rgt2[posu] :- rgt2[posb]) ) / (rgt2[rows(y)] - rgt2[1]) / 2
	cr = Ct1 + Ct2
	return(cr)
}
/* Computation of the Baseline Criterium for R2 */
real scalar quantilef::critbisf0()
{
	/* Get Gamma */
	gamm = 0

	/* common support condition */
	umin1 = max( (0, tg_1(0, gamm)) )
	umax1 = min( (1, tg_1(1, gamm)) )
	umin2 = max( (0, tg(0, gamm)) )
	umax2 = min( (1, tg(1, gamm)) )

	/* rank transformation */
	n1t = trunc( rows(y1) :* (umax1 - umin1) )
	n2t = trunc( rows(y2) :* (umax2 - umin2) )

	/* computation of threshold ranks */
	rangn1 = umin1 :+ ((range(0, nrange, 1)):/nrange) :* (umax1 :- umin1)
	tgrangn1 = tg( rangn1, gamm )
	rangn2 = umin2 :+ ((range(0, nrange, 1)):/nrange) :* (umax2 :- umin2)
	tg_rangn2 = tg_1( rangn2, gamm )
	w1 = weight1( rangn1 \ tg_rangn2 )
	w2 = weight2( tgrangn1 \ rangn2 )
	yn1 = quant( y1, rangn1 )
	qy2 = quant( y2, tgrangn1 )
	yn2 = quant( y2, rangn2 )
	qy1 = quant( y1, tg_rangn2 )
	y = ((yn1, qy2, rangn1, tgrangn1) \ (qy1, yn2, tg_rangn2, rangn2)), /*
		*/ ( w1, w2, (J(rows(yn1), 1 , cref1) \ J(rows(yn2), 1, cref2)) )

	/* Calculate Criteria */
	Cr = critb0()

	/* Criterium */
	return(Cr)
}
/*
** ==========================================================
** Misc
** ==========================================================
*/
/* Non-Parametric Computation of Quantile Derivatives, used for the test statistics */
real matrix quantilef::BlGart(real vector p, real vector yv)
{
	real matrix lambdp
	np = rows(p)
	n = rows(yv)
	m = trunc(sqrt(n))
	vk1 = (trunc(n :* p) :+ m) :* ((n :* p) :<= J(np, 1, n - m)) :+ n :* ((n :* p) :> J(np, 1, n - m))
	vk0 = (trunc(n :* p) :- m) :* ((n :* p) :>= J(np, 1, 1 + m)) :+ 1 :* ((n :* p) :< J(np, 1, 1 + m))
	lambdp = ( (yv[vk1] :- yv[vk0]) ) :/  ( n :* (vk1 :- vk0) )
	return(lambdp)
}
/* Computation of the Test Statistics */
real scalar quantilef::calctest(beta, gamm)
{
	real scalar stat
	umin1 = max( (0, tg_1(0, gamm)) )
	umax1 = min( (1, tg_1(1, gamm)) )
	umin2 = max( (0, tg(0, gamm)) )
	umax2 = min( (1, tg(1, gamm)) )
	rang1 = range(1, rows(y1), 1) :/ rows(y1)
	rang2 = range(1, rows(y2), 1) :/ rows(y2)
	tg_rang2 = tg_1(rang2, gamm)
	tgrang1 = tg(rang1, gamm)
	l1t = ( range(1, rows(y1), 1), (umin1 :<= rang1) :* (rang1 :<= umax1) )
	l2t = ( range(1, rows(y2), 1), (umin2 :<= rang2) :* (rang2 :<= umax2) )
	l1 = select( l1t, l1t[., 2] :!= 0 )
	l2 = select( l2t, l2t[., 2] :!= 0 )
	l1 = l1[., 1]
	l2 = l2[., 1]
	yn1 = y1[l1]
	yn2 = y2[l2]
	n1t = rows(yn1)
	n2t = rows(yn2)

	rangn1 = rang1[l1]
	tg_rangn2 = tg_rang2[l2]
	tgrangn1 = tgrang1[l1]
	rangn2 = rang2[l2]
	qy1 = quant(y1, tg_rangn2)
	qy2 = quant(y2, tgrangn1)
	y = ( (yn1, qy2, rangn1, tgrangn1) \ (qy1, yn2, tg_rangn2, rangn2) )
	((yn1, qy2) \ (qy1, yn2))
	(rangn1, tgrangn1) \ (tg_rangn2, rangn2)

	mu1 = 0.5 :* ( y[., 2] :+ zb(y[, 1], beta) )
	mu2 = 0.5 :* ( zb_1(y[., 2], beta) :+ y[., 1] )

	dlambd1 = BlGart( y[., 3], y1 )
	dlambd2 = BlGart( y[., 4], y2 )
	dmu1 = 0.5 :* ( dzb_u(y[1, 1], beta) :* dlambd1 :* (umax1 - umin1) :+ dtg_u(y[1, 4], gamm) :* dlambd2 :* (umax2 - umin2) )

	eps = 1e-8
	correc = ((umax2 - umin2) / (umax1 - umin1)) / ( (dzb_u(y[1, 2], beta)^2) / (dtg_u(y[1, 4], gamm) + eps) )
	vw1 = (1 :/ (dmu1:^2)) :/ (1 + correc)
	vw2 = vw1

	y = ( y, vw1, vw2, (J(rows(yn1), 1, cref1) \ J(rows(yn1), 1, cref2)) )

	cr = critb0(beta)
	stat = cr / ( (1 / rows(yn1)) + (1 / rows(yn2)) )
	return(stat)
}
end
/*---------- END Class quantilef ----------*/


/*
** ==========================================================
** Wrap-Arround Function for Optimization within Class
** ==========================================================
*/
version 12
mata:
/* wrt Beta */
void eval_critb(real scalar todo, real rowvector beta, class quantilef Q, real scalar v, real rowvector g, real matrix H)
{
	Q.critb(todo, beta, v, g, H)
}
/* wrt Gamma */
void eval_critbisf(real scalar todo, real rowvector beta, class quantilef Q, real scalar v, real rowvector g, real matrix H)
{
	Q.critbisf(todo, beta, Q, v, g, H)
}
end


/*
** ==========================================================
** STATA
** ==========================================================
*/
/*---------- START Command estquant ----------*/
*! version 1.20  11 July 2022
capture program drop estquant
program define estquant, sortpreserve eclass
	version 12
	syntax varlist [if] [in], /*
		*/ CATegory(varname) /*
		*/ [ /*
		*/ TRuncation /*
		*/ SHift /*
		*/ DIlation /*
		*/ INITR(real 0) /*
		*/ QRANGE(real 1000) /*
		*/ BVARiable(string) /*
		*/ BREPlication(real 50) /*
		*/ BSAMpling(real 100) /*
		*/ STRATA /*
		*/ OPTECHnique(string) /*
		*/ OPNMSimplexdeltas(real 1e-2) /*
		*/ MAXITeration(real 1e+3) /*
		*/ EPS1(real 1e-6) /*
		*/ EPS2(real 1e-6) /*
		*/ CI(string) /*
		*/ Level(real 95) /*
		*/ ]

	/* Variables */
	local vY `varlist'
	marksample touse
	markout `touse' `varlist' `category'
	tempname b V

	/* Check and Define Category */
	qui: tab `category'
	if( r(r) == 2 ){
		qui: sum `category'
		scalar cat1 = r(min)
		scalar cat2 = r(max)
	}
	else if( r(r) != 2 ){
		display as error "More than two categories are specified."
		exit 198
	}

	/* Check Variable */
	if( strpos("`vY'", " ") > 0 ){
		display as error "Multiple variables are not allowed."
		exit 103
	}

	/* Check Missing Values of Variable */
	qui: count if `vY' != . & `category' == cat1
	local cnt1_missing = r(N)
	qui: count if `vY' != . & `category' == cat2
	local cnt2_missing = r(N)
	if( `cnt1_missing' == 0 | `cnt2_missing' == 0 ){
		display as error "`vY' has no values."
		exit 198
	}

	/* Option: Truncation (Default Off) */
	local trnc = 0
	if( "`truncation'" != "" ){
		local trnc = 1
	}

	/* Option: Shift (Default Off) */
	local shft = 0
	if( "`shift'" != "" ){
		local shft = 1
	}

	/* Option: Dilation (Default Off) */
	local dltn = 0
	if( "`dilation'" != "" ){
		local dltn = 1
	}

	/* Option: Specify Initial Value of Truncation (Default Off) */
	local initrnc = 0
	scalar inivaltr = `initr'
	if( `initr' != 0 ){
		if( `initr' >= 0.5 ){
			display as error "Initial value of S in initr() must be lower than 0.5."
			exit 198
		}
		local initrnc = 1
	}

	/* Option: Confidence Interval */
	if( "`ci'" == "" ){
		/* Default: Normal-Based */
		local confi = "normal"
	}
	else if( "`ci'" != "" ){
		local confi = "`ci'"
	}
	if( "`confi'" != "bootstrap" & "`confi'" != "normal" ){
		display as error "Either bootstrap or normal must be specified in ci()."
		exit 198
	}

	/* Option: Confidence Interval Level */
	scalar cilevel = `level'
	if( ceil( (cilevel*10 - trunc(cilevel*10)) ) > 0 ){
		display as error "level() can have at most one digit after the decimal point"
		exit 198
	}

	/* Range */
	scalar nrange = `qrange'

	/* Bootstrap Option: strata */
	local bstrata = 0
	if( "`strata'" != "" ){
		/*bsample with strata option*/
		local bstrata = 1
	}

	/* Bootstrap Option: variable */
	if( "`bvariable'" == "" ){
		local bvar = "off"
	}
	else if( "`bvariable'" != "" ){
		local bvar = "`bvariable'"
	}

	/* Bootstrap Option: Error Check */
	if( "`bvar'" != "on" & "`bvar'" != "off" ){
		display as error "bvariable() requires either on or off."
		exit 198
	}

	/* Bootstrap Option: Error Check */
	if( `bstrata' == 1 & "`bvar'" == "on" ){
		display as text "note: strata option is ignored in the case of bvariable(on)."
	}

	/* Bootstrap Option: Replication and Sampling */
	local nbrep = `breplication'
	scalar nbsample = `bsampling'

	/* Bootstrap Option: Error Check for bvar option */
	if( "`bvar'" == "on" ){
		forvalues i = 1(1)`nbrep' {
			capture confirm new variable `vY'`i', exact
			if( _rc == 0 ){
				display as error "bvariable() requires variables named in a sequential order."
				exit 198
			}
			capture confirm new variable `category'`i', exact
			if( _rc == 0 ){
				display as error "bvariable() requires variables named in a sequential order."
				exit 198
			}
		}
	}

	/* Optimization Technique */
	if( "`optechnique'" == "" ){
		local optech = "nr"
	}
	else if( "`optechnique'" == "nr" | "`optechnique'" == "dfp" | "`optechnique'" == "bfgs" | "`optechnique'" == "nm" ){
		local optech = "`optechnique'"
		if( "`optechnique'" == "nm" ){
			display as text "Value of delta to built the simplex required by the Nelder–Mead method (`opnmsimplexdeltas')"
		}
	}
	else {
		display as error "optechnique() requires either nr, dfp, bfgs, or nm."
		exit 198
	}

	/* Optimization Option for the NM Method */
	scalar opnmsdeltas = `opnmsimplexdeltas'

	/* Convegence Criteria in Optimization */
	scalar maxIt = `maxiteration'
	scalar eps1 = `eps1'
	scalar eps2 = `eps2'

	/* Display Results of Each Iteration: 0: OFF, 1: ON: */
	local disp = 0

	/* Bootstrap Replications */
	if( `nbrep' == 0 ){
		display as text "Bootstrap replications skipped."
		/* Dummy Matrix for Bootstrap Replications */
		scalar nbsample = .
		matrix mP = J(1, 3, 0)
		matrix colnames mP = A D S
	}
	else if( `nbrep' == 1 ){
		display as error "Bootstrap requires 2 replications or more."
		exit 198
	}
	else if( `nbrep' >= 2 ){
		matrix mP = J(`nbrep', 3, 0)
		matrix colnames mP = A D S
		display as text "Bootstrap replications (`nbrep')
		if( `trnc' == 1 | `shft' == 1 | `dltn' == 1 ){
			/* Bootstrap Case 1: bvariable option off */
			if( "`bvar'" == "off" ){
				forvalues i = 1(1)`nbrep' {
					preserve
					local bcnt = 0
					local bconv = 0
					/* BOOTSTRAP WITH ERROR CHECK */
					while `bconv' == 0 {
						if( `bstrata' == 0 ){
							bsample round((nbsample/100) * _N)
						}
						else if( `bstrata' == 1 ){
							bsample round((nbsample/100) * _N), strata(`category')
						}
						/* ESTIMATION */
						mata: estquant("`vY'", "`category'", `trnc', `shft', `dltn', `initrnc', "`touse'", "`b'", "`V'", st_matrix("mP"), "`confi'", `disp')
						/* scalar: 0: not converged, 1: converged */
						local bconv = converged
						local bcnt = `bcnt' + 1
						if( `bcnt' > 50 ){
							/* If convergence not achieved */
							disp as error "{c |}Convergence Err{c |}"
							matrix vp = J(1, 3, .)
						}
					}
					matrix vP = vp
					matrix mP[`i', 1] = vP[1, 1]
					matrix mP[`i', 2] = vP[1, 2]
					matrix mP[`i', 3] = vP[1, 3]
					restore

					/* Display Iteration Process */
					if( `i' == 1 ){
						disp as text "{c TT}{hline 15}{c TT}"
					}
					if( `i' == trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  10%{c |}"
					}
					else if( `i' == 2*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  20%{c |}"
					}
					else if( `i' == 3*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  30%{c |}"
					}
					else if( `i' == 4*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  40%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) ){
						disp as text "{c |}Completed:  50%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  60%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 2*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  70%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 3*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  80%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 4*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  90%{c |}"
					}
					else if( `i' == `nbrep' ){
						disp as text "{c |}Completed: 100%{c |}"
						disp as text "{c BT}{hline 15}{c BT}"
					}
				}
			}
			/* Bootstrap Case 2: bvariable option on */
			if( "`bvar'" == "on" ){
				forvalues i = 1(1)`nbrep' {
					scalar nbsample = .
					/* ESTIMATION FOR BOOTSTRAP */
					mata: estquant("`vY'`i'", "`category'`i'", `trnc', `shft', `dltn', `initrnc', "`touse'", "`b'", "`V'", st_matrix("mP"), "`confi'", `disp')
					local bconv = converged
					if( `bconv' != 1 ){
						/* If convergence not achieved */
						disp as error "{c |}Convergence Err{c |}"
						matrix vp = J(1, 3, .)
					}
					/* STORE */
					matrix vP = vp
					matrix mP[`i', 1] = vP[1, 1]
					matrix mP[`i', 2] = vP[1, 2]
					matrix mP[`i', 3] = vP[1, 3]

					/* Display Iteration Process */
					if( `i' == 1 ){
						disp as text "{c TT}{hline 15}{c TT}"
					}
					if( `i' == trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  10%{c |}"
					}
					else if( `i' == 2*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  20%{c |}"
					}
					else if( `i' == 3*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  30%{c |}"
					}
					else if( `i' == 4*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  40%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) ){
						disp as text "{c |}Completed:  50%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  60%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 2*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  70%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 3*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  80%{c |}"
					}
					else if( `i' == trunc(`nbrep'/2) + 4*trunc(`nbrep'/10) ){
						disp as text "{c |}Completed:  90%{c |}"
					}
					else if( `i' == `nbrep' ){
						disp as text "{c |}Completed: 100%{c |}"
						disp as text "{c BT}{hline 15}{c BT}"
					}
				}
			}
		}
	}


	/* Display Main Results: 0: OFF, 1: ON: */
	local disp = 1

	/* Mata Main Program for Estimation */
	mata: estquant("`vY'", "`category'", `trnc', `shft', `dltn', `initrnc', "`touse'", "`b'", "`V'", st_matrix("mP"), "`confi'", `disp')

	/* ereturn in Stata */
	matrix colnames `b' = "Shift A" "Dilation D" "Truncation S"
	matrix colnames `V' = "Shift A" "Dilation D" "Truncation S"
	matrix rownames `V' = "Shift A" "Dilation D" "Truncation S"
	ereturn clear
	ereturn post `b' `V', esample(`touse')
	ereturn matrix B mP
	ereturn scalar N = nobs
	ereturn scalar N1 = nobs_c1
	ereturn scalar N2 = nobs_c2
	ereturn scalar qrange = nrange
	ereturn scalar brep = nbrep
	ereturn scalar bsample = nbsample
	ereturn scalar r2 = r2
	ereturn scalar cri = cri
	ereturn scalar cr1 = cr1
	ereturn scalar cr2 = cr2
	ereturn scalar mean_y = meany
	ereturn scalar mean_y1 = meany1
	ereturn scalar mean_y2 = meany2
	ereturn scalar sd_y = sdy
	ereturn scalar sd_y1 = sdy1
	ereturn scalar sd_y2 = sdy2
	ereturn local varname "`vY'"
	ereturn local cmd "estquant"


	/* Drop Variables and Scalars */
	scalar drop cat1 cat2 nrange nobs nobs_c1 nobs_c2 meany meany1 meany2 sdy sdy1 sdy2
	scalar drop nbrep nbsample converged maxIt eps1 eps2 inivaltr
end


/* ==========================
** Mata Main Program
** Call quantdist Class
========================== */
version 12
mata:
void estquant(vY, category, trnc, shft, dltn, initrnc, touse, bname, Vname, mP, ci, disp)
{
	/* Load Class (Mata Library) */
	class quantilef scalar Q

	/* Get Values */
	cref1 = st_numscalar("cat1")
	cref2 = st_numscalar("cat2")
	inivaltr = st_numscalar("inivaltr")
	nrange = st_numscalar("nrange")
	nbrep = rows(mP)
	nbsample = st_numscalar("nbsample")
	cilevel = st_numscalar("cilevel")
	opnmsdeltas = st_numscalar("opnmsdeltas")
	maxIt = st_numscalar("maxIt")
	eps1 = st_numscalar("eps1")
	eps2 = st_numscalar("eps2")
	optechnique = st_local("optechnique")

	/* Make Two Variables for Each Category */
	svarlist = category + " " + vY
	st_view(mY, ., svarlist, touse)
	mY1 = select(mY, mY[., 1] :== cref1)
	mY2 = select(mY, mY[., 1] :== cref2)
	y1 = sort( mY1[., 2], 1 )
	y2 = sort( mY2[., 2], 1 )
	y1 = select( y1, y1 :!= . )
	y2 = select( y2, y2 :!= . )

	/* Set Parameter: beta = (D, A) */
	beta_ini = J(1, (dltn + shft), 0)
	if( dltn == 1 ){
		beta_ini[1] = 1
	}

	/* Initial Setting */
	Q.init_setting_range( nrange )
	Q.init_setting_data( y1, y2 )
	Q.init_setting_category( (cref1, cref2) )
	Q.init_setting_spec( (trnc, shft, dltn) )
	Q.init_setting_params( beta_ini )
	Q.init_setting_conv( maxIt, eps1, eps2 )

	/* Baseline Criteria: (S, D, A) = (0, 1, 0) */
	pgammf = Q.gammapp_1(0)
	Cr0 = Q.critbisf0()

	/* Estimation of Gamma and Beta */
	/* Case: S != 0 */
	if( Q.truncation == 1 ){
		/* Auto-Specified Initial Value */
		if( initrnc == 0 ){
			pgammi = -2.0
			pgamm0 = pgammi
			Q.critbisf(0, pgammi, Q, voirs, g, H)
			/* Grid Search */
			for(i = 1; i <= 10; i++){
				pgammi = pgammi + 0.5
				Q.critbisf(0, pgammi, Q, voiri, g, H)
				if( voiri < voirs ){
					pgamm0 = pgammi
					voirs = voiri
				}
			}
		}
		/* Pre-Specified Initial Value */
		else if( initrnc == 1 ){
			gammi = inivaltr / (inivaltr - 1)
			pgammi = Q.gammapp_1( gammi )
			pgamm0 = pgammi
			Q.critbisf(0, pgammi, Q, voirs, g, H)
		}

		/* Optimization wrt gamma and beta */
		/* QUESTION: Why each result is different in nr algorithm? */
		/* Find initial values for the NR method */
		/* Repeat this step if the result is unsteady */
		if( optechnique == "nr" ){
			Sg = optimize_init()
			optimize_init_which(Sg, "min")
			optimize_init_evaluator(Sg, &eval_critbisf())
			optimize_init_evaluatortype(Sg, "d0")
			optimize_init_technique(Sg, "nm")
			optimize_init_nmsimplexdeltas(Sg, 0.1)
			optimize_init_argument(Sg, 1, Q)
			optimize_init_params(Sg, pgamm0)
			optimize_init_tracelevel(Sg, "none")
			optimize_init_conv_ptol(Sg, 1e-20)
			optimize_init_conv_vtol(Sg, eps1)
			optimize_init_conv_nrtol(Sg, eps2)
			optimize_init_conv_maxiter(Sg, maxIt)
			optimize_init_conv_warning(Sg, "off")
			pgamm0 = optimize(Sg)
			voirs = optimize_result_value(Sg)
			bconv = optimize_result_converged(Sg)
			st_numscalar("converged", bconv)
		}
		/* Optimization wrt gamma and beta */
		{
			Sg = optimize_init()
			optimize_init_which(Sg, "min")
			optimize_init_evaluator(Sg, &eval_critbisf())
			optimize_init_evaluatortype(Sg, "d0")
			if( optechnique == "nr" | optechnique == "dfp" | optechnique == "bfgs" ){
				optimize_init_technique(Sg, optechnique)
				optimize_init_singularHmethod(Sg, "hybrid")
			}
			else if( optechnique == "nm" ) {
				optimize_init_technique(Sg, optechnique)
				optimize_init_nmsimplexdeltas(Sg, opnmsdeltas)
			}
			optimize_init_argument(Sg, 1, Q)
			optimize_init_params(Sg, pgamm0)
			optimize_init_tracelevel(Sg, "none")
			optimize_init_conv_ptol(Sg, 1e-20)
			optimize_init_conv_vtol(Sg, eps1)
			optimize_init_conv_nrtol(Sg, eps2)
			optimize_init_conv_maxiter(Sg, maxIt)
			optimize_init_conv_warning(Sg, "off")
			pgamm0 = optimize(Sg)
			voirs = optimize_result_value(Sg)
			bconv = optimize_result_converged(Sg)
			st_numscalar("converged", bconv)
		}
		pgammf = pgamm0

		/* Store beta and gamma*/
		if( Q.dilation == 1 | Q.shift == 1 ){
			betaf = Q.betaf
		}
		gammf = Q.gammapp(pgammf)
	}
	/* Case: S = 0 */
	else if( Q.truncation == 0 ){
		pgammf = Q.gammapp_1(0)

		/* Optimization wrt beta */
		Q.critbisf(0, pgammf, Q, Cr, g, H)
		st_numscalar("converged", 1)

		/* Store beta and gamma*/
		if( Q.dilation == 1 | Q.shift == 1 ){
			betaf = Q.betaf
		}
		gammf = 0
	}

	/* Store Criterium */
	Cr1 = Q.Ct1
	Cr2 = Q.Ct2
	Cri = Q.Cr

	/* Computation of the Test Statistics */
	/*
	stat1 = Q.calctest(betaf, gammf)
	stat1
	*/

	/* Parameters */
	beta = J(2, 1, 0)
	gamm_e = - gammf / (1 - gammf)

	/* Results Table */
	if( Q.shift == 1 & Q.dilation == 1 ){
		beta[1] = betaf[1]
		beta[2] = betaf[2]
	}
	else if( Q.shift == 0 & Q.dilation == 1 ){
		beta[1] = betaf
		beta[2] = 0
	}
	else if( Q.shift == 1 & Q.dilation == 0 ){
		beta[1] = 1
		beta[2] = betaf
	}
	else {
		beta[1] = 1
		beta[2] = 0
	}
	vp = (beta[2], beta[1], gamm_e)
	sPName = ("Shift A", "Dilation D", "Truncation S")
	vA = sort(mP[.,1], 1)
	vD = sort(mP[.,2], 1)
	vS = sort(mP[.,3], 1)
	sdA = sqrt(variance(vA))
	sdD = sqrt(variance(vD))
	sdS = sqrt(variance(vS))
	zA = vp[1] / sdA
	zD = (vp[2] - 1) / sdD
	zS = vp[3] / sdS
	pA = 2 * (1 - normal(abs(zA)))
	pD = 2 * (1 - normal(abs(zD)))
	pS = 2 * (1 - normal(abs(zS)))
	alphaL = 0.5 * ((100 - cilevel) / 100)
	alphaU = 1 - 0.5 * ((100 - cilevel) / 100)
	lposu = ceil(nbrep * alphaL)
	lposb = floor(nbrep * alphaL)
	uposu = ceil(nbrep * alphaU)
	uposb = floor(nbrep * alphaU)
	if( lposb == 0 ){
		lposb = 1
	}
	if( ci == "bootstrap" ){
		cilA = 0.5 * (vA[lposu] + vA[lposb])
		cilD = 0.5 * (vD[lposu] + vD[lposb])
		cilS = 0.5 * (vS[lposu] + vS[lposb])
		ciuA = 0.5 * (vA[uposu] + vA[uposb])
		ciuD = 0.5 * (vD[uposu] + vD[uposb])
		ciuS = 0.5 * (vS[uposu] + vS[uposb])
	}
	else if( ci == "normal" ){
		cilA = vp[1] + sdA * invnormal(alphaL)
		cilD = vp[2] + sdD * invnormal(alphaL)
		cilS = vp[3] + sdS * invnormal(alphaL)
		ciuA = vp[1] + sdA * invnormal(alphaU)
		ciuD = vp[2] + sdD * invnormal(alphaU)
		ciuS = vp[3] + sdS * invnormal(alphaU)
	}
	b = vp
	V = diag((variance(vA), variance(vD), variance(vS)))
	mStat = (vp[1], sdA, zA, pA, cilA, ciuA ) \ /*
		*/ (vp[2], sdD, zD, pD, cilD, ciuD ) \ /*
		*/ (vp[3], sdS, zS, pS, cilS, ciuS )
	r2 = 1 - (Cri / Cr0)

	/* For Skipped Bootstrap Procedure */
	if( nbrep == 1 ){
		nbrep = 0
		V = J(3, 3, 0)
	}

	/* ereturn for Stata */
	st_matrix(bname, b)
	st_matrix(Vname, V)
	st_matrix("vp", vp)
	st_numscalar("nobs", Q.nobs)
	st_numscalar("nobs_c1", Q.nobs_c1)
	st_numscalar("nobs_c2", Q.nobs_c2)
	st_numscalar("meany", Q.mean_y)
	st_numscalar("meany1", Q.mean_y1)
	st_numscalar("meany2", Q.mean_y2)
	st_numscalar("sdy", Q.sd_y)
	st_numscalar("sdy1", Q.sd_y1)
	st_numscalar("sdy2", Q.sd_y2)
	st_numscalar("r2", r2)
	st_numscalar("nrange", nrange)
	st_numscalar("nbrep", nbrep)
	st_numscalar("nbsample", nbsample)
	st_numscalar("cr1", Cr1)
	st_numscalar("cr2", Cr2)
	st_numscalar("cri", Cri)

	/* Variable Name: Category Variable */
	scategory = abbrev(category, 16)

	/* Variable Name: Distribution Variable */
	sY = abbrev(vY, 16)

	/* Display Results */
	if( disp == 1 ){
		/* Summary Statistics */
		printf("\n\n")
		printf("{txt}%16s {c |}{space 5}Obs. {space 6}Mean {space 6}S.D. {space 7} Number of obs = {res}%8.0f \n", scategory, Q.nobs)
		printf("{hline 17}{c +}{hline 31} {txt}{space 7} Replications  = {res}%8.0f \n", nbrep)
		printf("{txt}%16.0f {c |}{res} %8.0f %11.4g %9.4g {space 8}{txt}BS Sample (%%) = {res}%8.4g \n", cref1, Q.nobs_c1, Q.mean_y1, Q.sd_y1, nbsample)
		printf("{txt}%16.0f {c |}{res} %8.0f %11.4g %9.4g {space 8}{txt}Pseudo R2     = {res}%8.4f \n", cref2, Q.nobs_c2, Q.mean_y2, Q.sd_y2, r2)
		printf("{hline 17}{c +}{hline 31} \n")
		printf("{txt}{space 11}Total {c |}{res} %8.0f %11.4g %9.4g \n", Q.nobs, Q.mean_y, Q.sd_y)
		printf("\n")
		printf("{hline 17}{c TT}{hline 64} \n")
		if( ci == "bootstrap" ){
			printf("{txt}{space 16} {c |}{space 3}Observed{space 3}Bootstrap{space 23}Bootstrap-based \n")
		}
		else if( ci == "normal" ){
			printf("{txt}{space 16} {c |}{space 3}Observed{space 3}Bootstrap{space 25}Normal-based \n")
		}
		printf("{txt}%16s {c |}{space 6}Coef.{space 3}Std. Err.{space 6}z{space 4}P>|z|{space 3}[%3.1f%% Conf. Interval]\n", sY, cilevel)
		printf("{hline 17}{c +}{hline 64} \n")
		if( Q.shift == 1 & Q.dilation == 1 & Q.truncation == 1 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[1], mStat[1,1], mStat[1,2], mStat[1,3], mStat[1,4], mStat[1,5], mStat[1,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[2], mStat[2,1], mStat[2,2], mStat[2,3], mStat[2,4], mStat[2,5], mStat[2,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[3], mStat[3,1], mStat[3,2], mStat[3,3], mStat[3,4], mStat[3,5], mStat[3,6] )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}nothing \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 1 & Q.dilation == 1 & Q.truncation == 0 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[1], mStat[1,1], mStat[1,2], mStat[1,3], mStat[1,4], mStat[1,5], mStat[1,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[2], mStat[2,1], mStat[2,2], mStat[2,3], mStat[2,4], mStat[2,5], mStat[2,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.6g %11.6g\n",
				sPName[3], mStat[3,1],.,.,.,.,. )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}S = 0 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 0 & Q.dilation == 1 & Q.truncation == 1 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1], mStat[2,2], mStat[2,3], mStat[2,4], mStat[2,5], mStat[2,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3], mStat[3,1], mStat[3,2], mStat[3,3], mStat[3,4], mStat[3,5], mStat[3,6] )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}A = 0 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 1 & Q.dilation == 0 & Q.truncation == 1 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1], mStat[1,2], mStat[1,3], mStat[1,4], mStat[1,5], mStat[1,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3], mStat[3,1], mStat[3,2], mStat[3,3], mStat[3,4], mStat[3,5], mStat[3,6] )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}D = 1 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 1 & Q.dilation == 0 & Q.truncation == 0 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1], mStat[1,2], mStat[1,3], mStat[1,4], mStat[1,5], mStat[1,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3], mStat[3,1],.,.,.,.,. )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}D = 1, S = 0 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 0 & Q.dilation == 1 & Q.truncation == 0 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1], mStat[2,2], mStat[2,3], mStat[2,4], mStat[2,5], mStat[2,6] )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3],mStat[3,1],.,.,.,.,. )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}A = 0, S = 0 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 0 & Q.dilation == 0 & Q.truncation == 1 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3], mStat[3,1], mStat[3,2], mStat[3,3], mStat[3,4], mStat[3,5], mStat[3,6] )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}A = 0, D = 1 \n")
			printf("{hline 82} \n")
		}
		else if( Q.shift == 0 & Q.dilation == 0 & Q.truncation == 0 ){
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[1], mStat[1,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[2], mStat[2,1],.,.,.,.,. )
			printf("{txt}%16s {c |}{res}%11.7g %10.5g   %6.2f %7.3f %12.7g %11.7g\n",
				sPName[3], mStat[3,1],.,.,.,.,. )
			printf("{hline 17}{c BT}{hline 64} \n")
			printf("{txt}Constrained parameters: {res}A = 0, D = 1, S = 0 \n")
			printf("{hline 82} \n")
		}
	}
}
end
/*---------- END Command estquant ----------*/