*! version 1.2.1 PR 01feb2008
program define fracdydx, rclass
	version 8
	local fp `e(fp_cmd2)'
	syntax [varlist(default=none max=1)] [ , Gen(str) GENSe(string) Powers(string) Coeffs(string) Deriv(int 1) ]
	if `deriv'<1 {
		di in red "invalid `deriv'"
		exit 198
	}
	tempvar x
	if "`fp'"=="fracpoly" {
		if "`e(fp_xpx)'"!="" {
			di in red "exponential transformation not supported"
			exit 198
		}
		if "`powers'`coeffs'"=="" {
			if "`varlist'"!="" {
				qui gen double `x'=`varlist' if e(sample)
			}
/*
			else qui gen double `x'=(`e(fp_x1)'+ /*
			 */ e(fp_shft))/e(fp_sfac) if e(sample)
*/
			else qui gen double `x'=(`e(fp_x1)'+e(fp_shft))/e(fp_sfac)
		}
		else local fp
	}
/*
	Store powers in macros and coefficients in locals
*/
	if "`fp'"!="fracpoly" {
		if "`varlist'"=="" {
			di in red "varname required"
			exit 198
		}
		if "`powers'"=="" {
			di in red "powers() required"
			exit 198
		}
		if "`coeffs'"=="" {
			di in red "coeffs() required"
			exit 198
		}
		cap local np=colsof(`powers')
		local rc=_rc
		if `rc' {
			// Treat as a numlist
			cap numlist "`powers'"
			if _rc {
				noi di as err `"invalid powers(`powers')"'
				exit 198
			}
			local powers `r(numlist)'
			local np: word count `powers'
		}
		forvalues i=1/`np' {
			if `rc' {
				local p`i': word `i' of `powers'
			}
			else local p`i'=`powers'[1,`i']
		}
		cap local r=colsof(`coeffs')
		local rc=_rc
		if `rc' {
			local r: word count `coeffs'
		}
		if `r'!=`np' {
			di in red "numbers of powers and coefficients differ"
			exit 198
		}
		forvalues i=1/`np' {
			if `rc' {
				local b`i': word `i' of `coeffs'
			}
			else local b`i'=`coeffs'[1,`i']
		}
		qui gen double `x'=`varlist'
	}
	else {
		local pwrs `e(fp_k1)'
		local np: word count `pwrs'
		forvalues i=1/`np' {
			local p`i': word `i' of `pwrs'
			local X`i': word `i' of `e(fp_xp)'
			local b`i'=_b[`X`i'']
		}
	}
	tempname small
	scalar `small'=1e-6
/*
	Compute powers and coefficients for (successive) derivative(s)
*/
	local m `np'
	local np0 `np'
	forvalues d=1/`deriv' {
		local i=`np'+1
		local p`i' .
		local c0 0
		local j1 0
		local j 1
		local pi .
		while `j'<=`np' {
			local pj `p`j''
			local jplus1=`j'+1
			local isrepeat=(abs(`p`jplus1''-`pj')<`small')	/* next power is repeated power */
			local islog=(abs(`pj')<`small')			/* this power is log */
			if abs(`pj'-`pi')>`small' {			/* pj is new power (not a repeated power) */
				local j0 `j'
			}
			if `islog' {
				local k=`j'-`j0'+1			/* deals with repeated powers of 0 */
			}
			else local k `pj'				/* power term */
			local c=`k'*`b`j''
			if `isrepeat' & !`islog' {			/* repeated power not of 0 */
				local c=`c'+(`j'-`j0'+1)*`b`jplus1''
			}
			if abs(`pj'-1)>`small'{
				local jj1=`j'-`j1'
				local c`jj1' `c'
				local p`jj1'=`pj'-1
			}
			else {
				if `j'==`j0' {
					local c0 `c'
					local m=`m'-1
					local j1 1
				}
				else {
					local jm1=`j'-1
					local c`jm1' `c'
					local p`jm1' 0
				}
			}
			local pi `pj'
			local j `jplus1'
		}
		forvalues j=1/`m' {
			local b`j' `c`j''
		}
		local np `m'
	}
/*
	Calc derivative
*/
	quietly {
		if "`gen'"!="" {
			tempvar new
			qui gen double `new'=`c0' if `x'!=.
		}
		tempname s4
		if `np'>0 {
/*
	Generate fractional powers.
*/
			tempvar lnx
			gen double `lnx'=log(`x')
			matrix `s4'=J(1,`np',0)
			local s2 `p1'
			local plast 0
			local hlast 1
			forvalues i=1/`np' {
				tempvar h`i'
				local pi `p`i''
				if abs(`pi'-`plast')>`small' {
					if abs(`pi')<`small' {
						gen double `h`i''=`lnx'
					}
					else if abs(`pi'-1)<`small' {
						gen double `h`i''=`x'
					}
					else gen double `h`i''= cond(`x'==0,0,`x'^`pi')
				}
				else {
					gen double `h`i''=`lnx'*`hlast'
				}
				local hlast `h`i''
				local plast `pi'
				if "`gen'"!="" {
					replace `new'=`new'+`c`i''*`h`i''
				}
				matrix `s4'[1,`i']=`c`i''
				if `i'>1 {
					local s2 `s2' `pi'
				}
			}
		}
	}
	if `np'<`np0' {
		local s2 `s2' .
	}
	if "`gen'"!="" {
		if "`varlist'"=="" {
			qui replace `new'=`new'/(e(fp_sfac)^`deriv')
		}
		cap drop `gen'
		rename `new' `gen'
		lab var `gen' "Derivative `deriv'"
	}

	if "`gense'"!="" & "`varlist'"=="" {
		// Calc SE of derivative
		tempname B VCE tmp
		local sfac = e(fp_sfac)
		matrix `B' = e(b)
		matrix `B' = `B'[1, "`X1'".."`X`np''"]	// submatrix for X's only
		matrix `VCE' = e(V)
		matrix `VCE' = `VCE'["`X1'".."`X`np''", "`X1'".."`X`np''"]	// submatrix for X's only
		// Replace x's temporarily with termise derivatives
		forvalues j=1/`np' {
			tempvar temp`j'
			rename `X`j'' `temp`j''
			qui gen double `X`j'' = `c`j''*`h`j''/_b[`X`j'']	// gives termwise derivative
		}
		_estimate hold `tmp'
		ereturn post `B' `VCE'
		cap drop `gense'
		_predict `gense' if !missing(`x'), stdp
		replace `gense' = `gense'/(`sfac'^`deriv')
		lab var `gense' "SE(derivative `deriv')"
		_estimate unhold `tmp'
		forvalues j=1/`np' {
			drop `X`j''
			rename `temp`j'' `X`j''
		}
	}
	return local powers `s2'
	return matrix coeffs `s4'
end