*! 1.0.0 NJC 27 February 2025
* lmoments 6.0.0 NJC 3 October 2012 
* 1.0.0 NJC 17 September 1997
* based on lshape v 2.0.1 PR 06Oct95.
program _lmoments, rclass    
	version 10 
	syntax varname(numeric) [if] [in] [, lmax(int 4) ]

	marksample touse

	forval l = 1/`lmax' { 
		local call `call' l_`l' 
	}

	mata: _lmo("`varlist'", "`touse'", "`call'", `lmax')  
  			
    return scalar N = scalar(N)  

	forval l = 1/`lmax' { 
		return scalar l_`l' = scalar(l_`l')
	}
     
	if `lmax' > 1 return scalar t = scalar(l_2) / scalar(l_1)
		
	forval l = 3/`lmax' { 
		return scalar t_`l' = scalar(l_`l') / scalar(l_2) 
	} 
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 lnames,  
real scalar lmax
) 
{ 
	real colvector x, result    
	real scalar n, j   

	x = st_data(., varname, usename) 
	n = length(x) 
	
	st_numscalar("N", n)  
	if (n == 0) return 

	_sort(x, 1)
	result = lmocoeff(n, lmax)' * x / n

	lnames = tokens(lnames) 
	for (j = 1; j <= length(lnames); j++) { 
		st_numscalar(lnames[j], result'[j])
	}
}

end