*! 6.0.0 NJC 3 October 2012
* 5.0.0 NJC 15 April 2010
* 4.0.1 NJC 16 October 2006
* 4.0.0 NJC 4 October 2006
* 3.3.0 NJC 20 April 2006
* 3.2.1 NJC 23 November 2004
* 3.2.0 NJC 18 October 2004
* 3.1.0 NJC 27 September 2004
* 3.0.0 NJC 1 July 2003
* 2.2.0 NJC 31 March 1999
* 2.1.0 NJC 12 February 1999
* 2.0.2 NJC 8 January 1999
* 2.0.1 NJC 1 September 1998
* 2.0.0 NJC 26 April 1998
* 1.0.0 NJC 17 September 1997
* based on lshape v 2.0.1 PR 06Oct95.
program lmoments, rclass byable(recall)
        version 10
        syntax [varlist] [if] [in] ///
        [, Format(str) lmax(numlist int >=4)  ALLobs variablenames short by(varlist) MISSing saving(str asis) * ]

	quietly {
		// screen out string variables
		ds `varlist', has(type numeric)
		local varlist "`r(varlist)'"

		// what to use
		if "`allobs'" != "" marksample touse, novarlist
		else marksample touse

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

		// variable(s) or group(s)
		local nvars : word count `varlist'
		local ng : word count `varlist'

		if "`by'" != "" {
			if `ng' > 1 {
				di as err ///
				"by() cannot be combined with `ng' variables"
				exit 198
			}

			tempvar group
			egen `group' = group(`by') if `touse', label `missing'
			su `group', meanonly
			local ng = r(max)
		}
		else {
			local group "`touse'"
			tokenize `varlist'
		}

		// initialisation
		if `ng' > _N {
			preserve
			set obs `ng'
			local preserved 1
		}
		else local preserved 0

		tempvar nnum nstr t myuse which
		tempname V

		capture label list which
		if _rc == 0 tempname mylbl
		else local mylbl "which"

		gen long `which' = _n
		compress `which'
		gen long `nnum' = .
		if "`lmax'" == "" local lmax = 4

		forval l = 1/`lmax' {
			tempvar l_`l'
			gen double `l_`l'' = .
                	label var `l_`l'' "l_`l'"
			local LTEMP `LTEMP' `l_`l''
		}

		gen double `t' = .
                label var `t' "t"

 		forval l = 3/`lmax' {
			tempvar t_`l'
			gen double `t_`l'' = .
                	label var `t_`l'' "t_`l'"
			local TTEMP `TTEMP' `t_`l''
		}

		gen byte `myuse' = `touse'

		// loop over variables(s) or groups(s)
		forval i = 1/`ng' {
			// get results
			if "`by'" != "" {
				replace `myuse' = `touse' & `group' == `i'
				mata: _lmo("`varlist'", "`myuse'", "`nnum'", "`LTEMP'", `i',`lmax')
			}
			else mata: _lmo("``i''", "`myuse'", "`nnum'", "`LTEMP'", `i', `lmax')

			// group or variable labels
			if "`by'" != "" {
				local V = trim(`"`: label (`group') `i''"')
			}
			else {
				local V = trim(`"`: variable label ``i'''"')
				if `"`V'"' == "" | "`variablenames'" != "" {
					local V "``i''"
				}
			}
			label def `mylbl' `i' `"`V'"', modify
		} // end loop over variables or groups

                // L-moment ratios
		replace `t'  = `l_2' / `l_1' in 1/`ng'
		forval l = 3/`lmax' {
			replace `t_`l'' = `l_`l'' / `l_2' in 1/`ng'
		}
	} // end quietly

	// -tabdisp- of results
        label val `which' `mylbl'
	if "`by'" != "" label var `which' "Group"
	else if "`allobs'" != "" label var `which' "Variable"
        else label var `which' "n = `rN'"

        if "`format'" == "" local format "%9.3f"
	if "`allobs'`by'" != "" {
		gen `nstr' = string(`nnum')
		label var `nstr' "n"
		local shown "`nstr'"
	}

	if "`short'" == "" {
	        tabdisp `which' in 1/`ng', ///
        	c(`shown' `l_1' `l_2' `l_3' `l_4') `options' f(`format')
	        tabdisp `which' in 1/`ng', ///
        	c(`t' `t_3' `t_4') `options' f(`format')
	}
        else {
                tabdisp `which' in 1/`ng', ///
                c(`shown' `l_1' `l_2' `t_3' `t_4') `options' f(`format')
        }

	// returned values
        ret scalar N = `nnum'[`ng']

	forval l = 1/`lmax' {
	        ret scalar l_`l' = `l_`l''[`ng']
	}

        ret scalar t = `t'[`ng']

	forval l = 3/`lmax' {
	        ret scalar t_`l' = `t_`l''[`ng']
	}

	// saving dataset?
	qui if `"`saving'"' != "" {
		if !`preserved' preserve
		keep in 1/`ng'
		keep `which' `nnum' `LTEMP' `t' `TTEMP'
		rename `which' which
		rename `nnum' n

		foreach v in `LTEMP' `t' `TTEMP' {
			rename `v' `: var label `v''
		}

		save `saving'
	}
end

mata :

real matrix bweights (real scalar n, real scalar k) {
	return(editmissing(comb((0::n-1), (0..k-1)) :/ comb(n-1, (0..k-1)), 0))
}

real matrix pweights(real scalar k) {
	real matrix w
	real scalar i, j
	w = J(k, k, .)

	for(i = 0; i < k; i++) {
		for(j = 0; j < k; j++) {
			w[i+1,j+1] = (-1)^(j-i) * exp(lnfactorial(j+i) - 2 * lnfactorial(i) - lnfactorial(j-i))
		}
	}

	return(editmissing(w, 0))
}

real matrix lmocoeff(real scalar n, real scalar k) {
	return(bweights(n, k) * pweights(k))
}

void _lmo(string scalar varname, string scalar usename, string scalar nname, string scalar lnames, real scalar i, real scalar lmax)
{
	real colvector x, result
	real scalar n

	x = st_data(., varname, usename)
	x = select(x, x :< .)
	n = length(x)
	st_store(i, nname, n)
	if (n == 0) return

	_sort(x, 1)
	result = lmocoeff(n, lmax)' * x / n
	if (n < lmax) result[(n + 1)::lmax] = J(lmax - n, 1, .)
	st_store(i, tokens(lnames), result')
}

end