*! lmose 1.0.0 NJC 13 April 2010
* lmoments 4.0.1 NJC 16 October 2006
* lmoments 1.0.0 NJC 17 September 1997
* lshape 2.0.1 PR 6 October 1995
program lmose, rclass byable(recall)
        version 10
        syntax varname(numeric) [if] [in] [, Format(passthru) * ]

	marksample touse
	qui count if `touse'
	if r(N) == 0 error 2000
 	if r(N) < 8 {
		di as txt "sample size too small for standard errors"
		exit 0
	}

	tempname V SE

	mata: _lmose("`varlist'", "`touse'", "`V'", "`SE'")

	if "`format'" == "" local format "format(%9.3f)"
	mat rownames `V' = l_1 l_2 l_3 l_4
	mat colnames `V' = l_1 l_2 l_3 l_4
	mat rownames `SE' = l_1 l_2 l_3 l_4 t t_3 t_4
	mat colnames `SE' = "  "

	di as txt _n "Variance matrix of sample L-moments" _c
	mat li `V', noheader `format' `options'
	di as txt _n "Standard errors of sample L-moments and ratios" _c
	mat li `SE', noheader `format' `options'

	ret matrix V = `V'
	ret matrix SE = `SE'
end

mata:

void _lmose(string scalar varname, string scalar usename,
string scalar Vname, string scalar SEname)
{
	real colvector x, i, b, work, IK, IL, JL, JK, prod, xixj, range, L, SE
	real scalar n, d, j, k, l, ikl, nn
	real matrix kl, V, C

	x = st_data(., varname, usename)
	x = select(x, x :< .)
	_sort(x, 1)
	n = length(x)
	i = (1..n)'
	b = J(4,1,0)
 	b[1] = mean(x)
	d = n
	work = x

	for (j = 1; j <= 3; j++) {
		d = d * (n - j)
		work = work :* (i :- j)
		b[j+1] = sum(work) / d
	}

	kl = J(0, 2, .)
	for (k = 0; k <= 3; k++) {
		for (l = k; l <= 3; l++) {
			kl = kl \ (k,l)
		}
	}

	V = J(4,4,0)

	for (ikl = 1; ikl <= rows(kl); ikl++) {
		k = kl[ikl, 1]
		l = kl[ikl, 2]

		JK = JL = IL = IK = J(n, 1, 1)
		for (j = 1; j <= k; j++) {
			IK = IK :* (i :- j)
			JK = JK :* (i :- (l+j+1))
		}

		for (j = 1; j <= l; j++) {
			IL = IL :* (i :- j)
			JL = JL :* (i :- (k+j+1))
		}

		nn = exp(lnfactorial(n) - lnfactorial(n-k-l-2))

		// want to work with pairs of order statistics x_(i) <= x_(j),
		// but we can avoid looping over both i and j

		prod = J(n, 1, 0)
		for (j = 1; j < n /* NB */; j++) {
			range = (j+1) \ n
			xixj = (x[j] / nn) * x[|range|]
			prod[|range|] = prod[|range|] +
			(IK[j] * JL[|range|] + IL[j] * JK[|range|]) :* xixj
		}
		V[k+1, l+1] = V[l+1, k+1] = sum(prod)
	}

	V = b * b' - V
	C = (1, 0, 0, 0\-1, 2, 0, 0\1, -6, 6, 0\-1, 12, -30, 20)
	L = C * b
 	V = C * V * C'

	SE = J(7,1,0)
	SE[|1\4|] = sqrt(diagonal(V[|(1,1)\(4,4)|]))
	SE[5] = V[2,2]/(L[2] * L[2]) + V[1,1]/(L[1] * L[1]) - 2 * V[1,2]/(L[1] * L[2])
	SE[5] = sqrt(SE[5] * (L[2]/L[1])^2)
	SE[6] = V[3,3]/(L[3] * L[3]) + V[2,2]/(L[2] * L[2]) - 2 * V[2,3]/(L[2] * L[3])
	SE[6] = sqrt(SE[6] * (L[3]/L[2])^2)
	SE[7] = V[4,4]/(L[4] * L[4]) + V[2,2]/(L[2] * L[2]) - 2 * V[2,4]/(L[2] * L[4])
	SE[7] = sqrt(SE[7] * (L[4]/L[2])^2)

	st_matrix(Vname, V)
	st_matrix(SEname, SE)
}

end