*! version 3.1.1 JPR 12-Dec-94.         (sg26.3: STB-25)
program define fpx
	version 3.1
	loc varlist "req ex min(2) max(3)"
	loc if "opt"
	loc in "opt"
	loc weight "aweight fweight"
	#delimit ;
	loc options "POwers(str) ADDpowers(str) FIXpowers(str) BAse(str)
CMd(str) CONTinuous(str) EXPx(str) NAme(str) DF(int 4) LOg
ORIgin(str) DEVthr(str) noCONST ZERo *" ;
	#delimit cr
	loc small 1e-6
	parse "`*'"
	if "`const'"!="" { loc options "`options' nocons" }
	if "`devthr'"!="" { conf num `devthr' }
	if "`powers'"=="1" { local df 1 }
	_jprfpin `df' "`cmd'" -1 "" /* gplot and saving not used */
	loc df $S_1
	loc cmd "$S_5"
	loc anova $S_6
	loc normal $S_7
	loc qreg $S_8
	loc glm $S_9
	loc blogit $S_10
	loc dist $S_11
	loc odddf=(2*int(`df'/2)<`df')
	if `df'==1 { local powers 1 }
	parse "`varlist'", parse(" ")
	gl S_E_cmd
	loc y `1'
	if `blogit' {
		loc cmd=substr("`cmd'", 2, .)
		loc mhs `2'
		loc rhs `3'
	}
	else {
		if "`3'"!="" { error 103 }
		loc rhs `2'
	}
	loc name=substr("`name'",1,6)
	if "`name'"=="`rhs'" {
		di in red "`name' already defined"
		exit 110
	}
	if `normal' {
	 loc devcal "\$_nobs*(1+log(2*_pi*_result(4)/\$_nobs)-\$_mnlnwt)"
	 loc pcal "fprob(\$_n1,\$_n2,(exp(\$_d/\$_nobs)-1)*(\$_n2)/(\$_n1))"
	}
	else {
		loc pcal "chiprob(\$_n1,\$_d)"
		if `dist'==2 {
			loc devcal "-2*\$S_E_ll"
		}
		else if `glm' {
			loc devcal "\$S_3"
		}
		else {
			loc devcal "-2*_result(2)"
		}
	}
/*
	Parse the powers
*/
	if "`powers'"=="" { loc powers "-2,-1,-.5,0,.5,1,2,3" }
	loc pwrs "`powers' `addpowe'"
	tempvar p
	qui gen `p'=.
	_jprfppp `p' "`pwrs'" 1 1
	loc np $S_1
	if `np'<1 {
		di in red "no powers given"
		exit 2002
	}
	loc pwrs "$S_3"
	tempvar tmp x lnx yl yh yfl yfh stdp
	qui {
		egen int `tmp'=rmiss(`y' `mhs' `base') `if' `in'
 		gen `x'=`rhs' if `tmp'==0
		drop `tmp'
	}
/*
	Deal with weights.
*/
	loc mnlnwt 0 /* mean log normalized weights */
	qui if "`exp'" != "" { 
		tempvar s
		gen `s' `exp'
		replace `x'=. if `s'<=0|`s'==.
		replace `s'=. if `x'==.
		if "`weight'"=="fweight" { 
			capture assert `s'==int(`s') if `x'!=.
			if _rc { noi error 401 } 
		}
		sum `s' 
		noi di "(sum of wgt is " _result(1)*_result(3) ")"
		replace `s'=log(`s'/_result(3))
		sum `s'
		loc mnlnwt=_result(3)
		drop `s'
		loc exp "[`weight'`exp']"
	}
	quietly {
		gen `lnx'=.
		loc lin=(`np'==1 & `p'[1]==1)
		if `lin' { loc zero }	
		_jprfpxo `x' `lnx' `lin' "`expx'" "`origin'" "`zero'"
		loc nobs $S_1
		loc zeta $S_2
		loc shift $S_3
		loc kx "$S_4"
		if "`devthr'"!="" {
			gen `yl'=.
			gen `yh'=.
			gen `yfl'=.
			gen `yfh'=.
			gen `stdp'=.
		}
	}
/*
	If `cmd' is blogit or bprobit,
	expand dataset to enable use of logit or probit command.
*/
	qui if `blogit' {
		_jprfplx `y' `mhs'
		loc y "_y_eq_1"
		loc exp "[freq=_pop]"
*		preserve
*		keep `y' `x' `mhs' `lnx' `p' `base'
*		drop if `x'==.
*		_jprfplx `y' `mhs'
*		drop `y'
*		rename _outcome `y'
*		loc exp "[freq=_pop]"
	}
/*
	Create variables corresponding to `fixpowers', if any.
*/
	if "`fixpowe'"!="" {
		tempvar xpp
		qui gen `xpp'=.
		_jprfppp `xpp' "`fixpowe'" 0 0
		if $S_2 {
			di in red "fixpowers must be in increasing order"
			exit 198
		}
		loc nxp $S_1
		loc j 1
		loc h
		while `j'<=`nxp' {
			loc hj "h`j'"
			tempvar `hj'
			qui gen ``hj''=.
			loc h "`h' ``hj''"
			loc j=`j'+1
		}
		_jprfpgn `x' `lnx' `xpp' `nxp' `h'
/*
	Determine positions of terms in `fixpowe'
*/
		loc j 1
		while `j'<=`np' {
			loc pj=`p'[`j']
			loc i 1
			loc pos 0
			while `i'<=`nxp' {
				if abs(`xpp'[`i']-`pj')<`small' {
					loc pos `i'
				}
				loc i=`i'+1
			}
			loc pos`j' `pos'
			loc j=`j'+1
		}
	}
	else {
		loc nxp 0
		loc j 1
		while `j'<=`np' {
			loc pos`j' 0
			loc j=`j'+1
		}
	}
/*
	Count residual and model df.
*/
	if `anova' {
		loc a "anova"
		if "`continu'"!="" { loc cts "cont(`continu')" }
	}
	else loc a "reg"
	qui `a' `y' `base' `exp' if `x'!=., `cts'
	loc rdf=_result(5)
	loc mdf=`nobs'-`rdf'
	if `rdf'<(1+`df'+`nxp') { error 2001 }
/*
	Calc deviance=-2(log likelihood) for base regression,
	allowing for possible weights.
*/
	if `anova' {
		if "`continu'"!="" { loc cts "cont(`continu')" }
	}
	qui `cmd' `y' `base' `exp' if `x'!=., `cts' `options'
	if `qreg' {
		loc dev0=`nobs'*log($S_E_msd/$S_E_rsd)
		loc rsd0 $S_E_rsd
	}
	else {
		loc dev0=`devcal'
		if `normal' { loc rsd0=_result(9) }
	}
	loc m=1+int((`df'-1)/2)
	if "`log'"!="" {
		di in gr "Dev. diff." _cont
		if `normal' { di in gr "  Res SD " _cont }
		loc i 1
		while `i'<=`m' {
			di in gr "   m = " `m'-`i'+1 _cont
			loc i=`i'+1
		}
		di _n
	}
/*
	Go!
*/
	loc i 1
	while `i'<=`np' {
		loc pi=`p'[`i']
		tempvar x`i'
		if `pos`i''>0 {
			qui gen `x`i''=`h`pos`i'''*`lnx'
		}
		else {
			qui gen `x`i''=cond(`pi'==0, `lnx', /*
			 */ cond(`pi'==1, `x', cond(`x'==0,0,`x'^`pi')))
		}
		loc i=`i'+1
	}
	loc i 0
	while `i'<`m' {
		loc j`i' 0
		loc i=`i'+1
		tempvar xp`i'
		qui gen `xp`i''=.
		loc j=`i'+`i'-1
		loc dev`j' `dev0' /* min deviance for df=j */
		loc j=`j'+1
		loc dev`j' `dev0'
	}
	loc j`m' 1
	loc i `m'
	loc kount 0
	loc deg `j`m''
	loc done 0
	while !`done' {
		if `i'<`m' {
			loc ji `j`i''
			qui replace `xp`i''=`x`ji''
			loc l `i'
			while `l'<`m' {
				loc l=`l'+1
				loc j`l' `ji'
				loc l1=`l'-1
				qui replace `xp`l''=`xp`l1''*`lnx'
			}
			if "`log'"=="" { di "." _cont }
		}
		else qui replace `xp`m''=`x`j`m'''
/*
	Test for any power being 1 (i.e. odd-df model)
*/
		loc one 0
		loc l 0
		while `l'<`m' {
			loc l=`l'+1
			if `p'[`j`l'']==1 { loc one 1 }
		}
* di "deg,df,jm=" `deg',`dfi',`j`m''
		if !`odddf' | (`deg'<`m') | `one' {
			loc dfi=2*`deg'-`one'
			loc i 1
			loc xlist
			while `i'<=`m' {
				if `j`i''>0 { loc xlist "`xlist' `xp`i''" }
				loc i=`i'+1
			}
			loc kount=`kount'+1
			if `anova' { loc cts "cont(`continu' `xlist' `h')" }
			qui `cmd' `y' `xlist' `h' `base' `exp', `cts' `options'
			if `qreg' {
				loc dev=`nobs'*log($S_E_msd/$S_E_rsd)
				loc rsd $S_E_rsd
			}
			else {
				loc dev=`devcal'
				if `normal' { loc rsd=_result(9) }
			}
			_jprfpdx "`devthr'" `dev' `yl' `yh' `yfl' `yfh' `stdp'
			if `dev'<`dev`dfi'' {
				loc dev`dfi' `dev'
				if `normal' { loc rsd`dfi' `rsd' }
				loc i 1
				loc pm`dfi'
				while `i'<=`m' {
					if `j`i''>0 {
						loc pm=`p'[`j`i'']
						loc pm`dfi' "`pm`dfi'' `pm'"
					}
					loc i=`i'+1
				}
			}
			if "`log'"!="" {
				di in ye %9.3f `dev0'-`dev' _cont
				if `normal' { di " " %8.0g `rsd' _cont }
				loc i 1
				while `i'<=`m' {
					if `j`i''==0 { di _skip(6) ". " _cont }
					else di in ye %8.1f `p'[`j`i''] _cont
					loc i=`i'+1
				}
				di
			}
		}
/*
	Increment the first possible index (of loop i) among indices of
	loops m, m-1, m-2, ..., 1
*/
		loc i `m'
/*
	Finish after all indices have achieved their upper limits (np).
*/
		while `j`i''==`np' {
			loc i=`i'-1
		}
		if `i'==0 { loc done 1 }
		else {
			if `j`i''==0 { loc deg = `m'-`i'+1 }
			loc j`i'=`j`i''+1
		}
	}
/*
	Update the results for even df to include odd df
*/
	loc i 2
	while `i'<=`df' {
		loc j=`i'-1
		if `dev`j''<`dev`i'' {
			loc dev`i' `dev`j''
			if `normal' { loc rsd`i' `rsd`j'' }
			loc pm`i' "`pm`j''"
		}
		loc i=`i'+2
	}

	if "`log'"=="" { di }
	loc dash "     --"
	loc d=64+15*`normal'
	di in gr _n in ye "`cmd'" in gr " models for " in ye "`y'" /*
		 */ in gr " on " in ye "`rhs'" /*
		 */ in gr " searched: " in ye `kount' in gr "."
	di in gr _n "Model" _col(14) "Deviance     " _cont
	if `normal' { di in gr "Residual SD   " _cont }
	di in gr " Akaike IC      Gain   ---- Powers ----" 
	di in gr _dup(`d') "-"
	di in gr "Base (*)" _col(13) in ye %9.3f `dev0' _cont
	if `normal' { di in ye _skip(5) %9.0g `rsd0' _cont }
	di in ye _skip(5) %9.3f `dev0'+2*`mdf' /*
		*/ in gr _skip(3) "`dash'`dash'"
	loc i 1
	while `i'<=`df' {
		loc mi=1+int((`i'-1)/2)
		di in gr "m = `mi' (`i' df)" _col(13) in ye %9.3f `dev`i'' _cont
		if `normal' { di in ye _skip(5) %9.0g `rsd`i'' _cont }
		di in ye _skip(5) %9.3f `dev`i''+2*(`mdf'+`i'+`nxp') /*
		 */ _skip(3) %9.3f `dev1'-`dev`i'' _skip(2) "`pm`i''"
		loc i=`i'+1
	}
	di in gr _dup(`d') "-"

	if "`base'"=="" { loc b "[none]" }
	else loc b "`base'"
	di _n in gr "(*) Base model = " in ye "`b'" /*
	 */	in gr " (" in ye `nobs' in gr " obs.)"
	if `nxp'>0 {
		di in gr "`plus' Fixed powers included in all FP models: " /*
		 */ in ye "`fixpowe'"
	}
	if "`zeta'"!="" {
		version 2.1
		local zeta = `zeta'	/* round */
		version 3.1
		if `shift' {
			if `zeta'<0 {
				local z "-"
				local zeta = -`zeta'
			}
			else local z "+"
			di in gr _n "X transformed to X`z'`zeta'."
			}
			else {
			if `zeta'>0 {
				local z "-"
			}
			else {
				local zeta = -`zeta'
				local z "+"
			}
		di in gr _n "X transformed to (X`z'`zeta')/(Xmax`z'`zeta')."
		}
	}
	if "`devthr'"!="" {
		cap drop _etal
		cap drop _etah
		cap drop _fl
		cap drop _fh
		rename `yl' _etal
		rename `yh' _etah
		rename `yfl' _fl
		rename `yfh' _fh
		rename `stdp' _stdp
	}
	if "`name'"=="" { loc name X }
	cap drop `name'
	rename `x' `name'
	lab var `name' "`rhs'"
	loc pwrs "`pm`df'' `fixpowe'"
	fpgen `name', powers(`pwrs') replace `zero'
	loc xp "$S_1"
	if `anova' { loc cts "cont(`continu' `xp')" }
	if "`log'"=="" { loc q "quietly" }
	`q' `cmd' `y' `xp' `base' `exp', `cts' `options'
	`q' di in gr "Fractional power(s) of `name' used: " in ye "`pwrs'"

	zap_s
	gl S_1 `nobs'
	gl S_2 `dev0'
	loc i 1
	while `i'<=`df' {
		loc j=`i'+2
		loc k=`j'+`df'
		gl S_`j' `dev`i''
		gl S_`k' "`pm`i''"
		loc i=`i'+1
	}
	gl S_23 `zeta'

	gl S_E_fprp "no"
	gl S_E_fp "fp"
	gl S_E_cmd "`cmd'"
	gl S_E_bloc `blogit'
	gl S_E_dist `dist'
	gl S_E_depv "`y'"
	gl S_E_x "`name'"
	gl S_E_xp "`xp'"
	gl S_E_xpx "`kx'"
	gl S_E_rhs "`rhs'"
	gl S_E_base "`base'"
	gl S_E_exp "`exp'"
	gl S_E_pwrs "`pwrs'"
	gl S_E_fixp "`fixpowe'"
	gl S_E_cmd "`cmd'"
end