/*******************************************************************************
Define Mata functions to be used in do-files as needed.
These functions are used for semiparametric estimations.
Updates
2013-05-28: allow multiple independent vars/weights in conditional expectation
calculation. Conditional expectation should be calculated for each
column of the y matrix and returned in a matrix of the same dimension.
*******************************************************************************/
*! lsls.do version 1.0 2014-10-11
*! author: Michael Barker mdb96@georgetown.edu
version 11
program drop _all
* Clear Mata
mata
mata clear
mata set matafavor speed
// mata set matalnum on
mata set matastrict on
/*******************************************************************************
Define CEXP parameter structure
*******************************************************************************/
struct cexp_parameters {
pointer(function) scalar method
pointer(function) scalar kernel
pointer(function) scalar dkernel
pointer(function) scalar bwidth
real scalar twice
real vector optbwidth
// Use these if it becomes necessary to include trim var in structure.
// real vector tx
// real scalar trim
}
/*******************************************************************************
Conditional Expectation Structure Initializer Fxns
*******************************************************************************/
// Define conditional expectation structure and save in referenced mopt. problem.
// Usage is specific to moptimize problems
// M: moptimize problem name.
// i: index for user info, where new cexp_parameters structure will be saved.
void mopt_cexp_define( transmorphic scalar M ,
real scalar i ,
pointer(function) vector fxnlist ) {
struct cexp_parameters scalar CE
CE = cexp_define(fxnlist)
moptimize_init_userinfo(M , i, CE)
}
// Define conditional expectation structure.
struct cexp_parameters scalar cexp_define(
| pointer(function) vector fxnlist
, real scalar twice ) {
struct cexp_parameters scalar CE
// Declare Semi-Parametric Parameters
// Default Parameters
if (args()==0) {
CE.method = &cexp_vec()
CE.kernel = &kernel_twicing()
CE.dkernel = &dkernel_twicing()
CE.bwidth = &bwidth_pilot()
CE.twice = 0
}
// User parameters
else {
CE.method = fxnlist[1]
CE.kernel = fxnlist[2]
CE.dkernel = fxnlist[3]
CE.bwidth = fxnlist[4]
CE.twice = twice
}
return(CE)
}
/*******************************************************************************
CEXP evaluator:
Pass conditional expectation input data and parameters to chosen evaluator
Return conditional expectation vector
Description of twicing as a discrete process from:
Smoothing Techniques for Curve Estimation
Lecture Notes in Mathematics Volume 757, 1979, pp 191-195
W Stuetzle, Y Mittal - Smoothing techniques for curve estimation, 1979
*******************************************************************************/
real matrix cexp( real matrix y, real matrix X, struct cexp_parameters CE,
| real vector tx, real vector h) {
// db "Conditional Expectation Called"
// If indicator trimming vector and bandwidth were not passed,
if (args()==3) {
// use a vector of ones. -> compute cond. exp. for all observations.
tx = J(rows(X),1,1)
// Calculate optimal bandwidth
h=(*CE.bwidth)(y,X,CE,tx)
}
// Trimming vector passed but bandwidth was not
else if (args()==4) {
// Calculate optimal bandwidth
h=(*CE.bwidth)(y,X,CE,tx)
}
// Update current bandwidth
CE.optbwidth=h
// calculate conditional expectation
real matrix cexp
cexp = (*CE.method)(y,X,CE,tx)
// Twicing: Discrete Process
if (CE.twice==1) {
real vector r, rhat
r = y-cexp
rhat = (*CE.method)(r,X,CE,tx)
cexp = cexp+rhat
}
return(cexp)
}
/*******************************************************************************
Kernel weighted distance using NxN matrix.
*******************************************************************************/
real matrix kdist(real matrix X , struct cexp_parameters scalar CE) {
real scalar n, k, j, h
real vector H, ones, z
real matrix K, Z
n = rows(X)
k = cols(X)
H = CE.optbwidth
// declare row vector of 1's to project the conditioning vector into nxn
ones = J(1,n,1)
// kernel weighted distance matrix
K = J(n,n,1)
// "leave one out" - Set diagonal of K to zero
_diag(K,0)
// loop through X's
for (j=1 ; j<=k ; j++) {
z = X[.,j]
h = H[j]
// Duplicate z into n columns
Z = z*ones
// Create square difference matrix with bandwidth factor
// update probability accumulation matrix
K = K :* ((*CE.kernel)((Z-Z'):/h) )
} // end for
// return matrix of kernel weighted distances
return(K)
}
/*******************************************************************************
Conditional Expectation: Uses NxN kernel distance weight matrix from kdist()
*******************************************************************************/
real matrix cexp_mat(real matrix y, real matrix X,
struct cexp_parameters scalar CE,
real vector tx) {
// Declarations
real matrix K
real matrix num
real vector den
// Get Kernel weighted distance matrix
K = kdist(X,CE)
// Calculate the conditional mean
num = K*y
den = quadrowsum(K)
return(tx :* (num :/ den))
}
/*******************************************************************************
Conditional Expectation Vector Method:
Compute kernel density by observation/vector
*******************************************************************************/
real matrix cexp_vec(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
real scalar i , k
real vector h
real matrix cexp , r
// Get bandwidth, row vector size K.
h = CE.optbwidth
// Initialize cexp accumulation vector:
// value=0 for obs that will be trimmed.
cexp = J(rows(X),cols(y),0)
// Calculate conditional expectation row by row.
for (i=1 ; i<=rows(X) ; i++) {
// Skip rows that will be trimmed.
if (tx[i]) {
// r Nxcols(x)
r = (*CE.kernel)((X[i,.] :- X):/h)
// Leave one out
r[i,1]=0
// Accumulate product across columns
// r[.,1] Nx1
for (k=2 ; k<=cols(X) ; k++) {
r[.,1] = r[.,1]:*r[.,k]
}
// cexp[i,.] 1xcols(y)
cexp[i,.] = cross(r[.,1],y) :/ quadsum(r[.,1])
}
}
return(cexp)
}
/*******************************************************************************
Derivative of Cond. Exp. WRT Y - Vector Method:
Compute derivative of kernel density by observation/vector
*******************************************************************************/
real matrix dcexpdy_vec(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
real scalar i , k
real vector h
real matrix cexp , r, dcexpdy
// Get bandwidth, row vector size K.
h = CE.optbwidth
// Initialize cexp accumulation vector:
// value=0 for obs that will be trimmed.
dcexpdy = J(rows(X),cols(y),0)
// Calculate conditional expectation row by row.
for (i=1 ; i<=rows(X) ; i++) {
// Skip rows that will be trimmed.
if (tx[i]) {
// r Nxcols(x)
r = (*CE.kernel)((X[i,.] :- X):/h)
// Leave one out
r[i,1]=0
// Accumulate product across columns
// r[.,1] Nx1
for (k=2 ; k<=cols(X) ; k++) {
r[.,1] = r[.,1]:*r[.,k]
}
// cexp[i,.] 1xcols(y)
dcexpdy[i,.] = cross(r[.,1],y) :/ quadsum(r[.,1])
}
}
return(dcexpdy)
}
/*******************************************************************************
Derivative of Conditional Expectation with respect to betas.
Betas are coefficients of conditional expectation Index: I = x1 + X2*beta
If additional indexes are included, pass them in optional argument, I2
Derivatives are taken only w.r.t. betas of first index.
Vector Method
*******************************************************************************/
real matrix dcexpdb_vec(real vector y , real colvector I1 , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx , real vector h ,| real matrix I2) {
real scalar i , k
real vector r1, dr1
real matrix dcexpdb , r2
// Get bandwidth, row vector size K.
// h = CE.optbwidth
// Adjust for derivative
// h = h :* rows(I1)^(2/99)
// h = h :^(9/11)
// Initialize cexp accumulation vector:
// value=0 for obs that will be trimmed.
dcexpdb = J(rows(X),cols(X),0)
// Initialize index
// I1 = X*b'
// Calculate conditional expectation row by row.
for (i=1 ; i<=rows(X) ; i++) {
// Skip rows that will be trimmed.
if (tx[i]) {
// r1 (x1i - x1j) Nx1
// subtract all rows of X from current row, col 1 only
r1 = (*CE.kernel)((I1[i] :- I1):/h[1])
dr1 = (*CE.dkernel)((I1[i] :- I1):/h[1]) :* (X[i,.] :- X) / h[1]
// Leave one out
r1[i] =0
dr1[i,.]=J(1,cols(X),0)
// calculate kernel weights of other indeces, if any
if (args()==7) {
// r (X2i - X2j) Nxcols(X)-1
r2 = (*CE.kernel)((I2[i,.] :- I2) :/ h[(2..cols(h))] )
// Accumulate product across columns
for (k=1 ; k<=(cols(r2)) ; k++) {
r1 = r1:*r2[.,k]
dr1 = dr1:*r2[.,k]
}
}
// cexp[i,.] 1xcols(y)
real scalar p, dp
real rowvector ydp, yp
p = quadsum(r1)
dp = quadcolsum(dr1)
yp = cross(y,r1)
ydp = cross(y,dr1)
dcexpdb[i,.] = (ydp:*p - yp:*dp) / p^2
}
}
return(dcexpdb)
}
/*******************************************************************************
Derivative of conditional expectation w.r.t. index
Change to dh format
*******************************************************************************/
real matrix dcexpdI_vec(real vector y , real matrix I1 ,
struct cexp_parameters scalar CE ,
real vector tx , real vector h ,
| real matrix I2) {
real scalar i , k
real vector r1, dr1, dcexpdI , p, yp, dLj, dKj
real matrix r2
// Get bandwidth, row vector size K.
// h = CE.optbwidth
// Adjust for derivative
// h = h :* rows(I1)^(2/99)
// h = h :^(9/11)
// Initialize cexp accumulation vector:
// value=0 for obs that will be trimmed.
dcexpdI = J(rows(y),1,0)
p = J(rows(y),1,0)
yp = J(rows(y),1,0)
dLj = J(rows(y),1,0)
// Calculate derivative of Li wrt Ii for each row
for (i=1 ; i<=rows(y) ; i++) {
// Skip rows that will be trimmed.
if (tx[i]) {
// r1 (x1i - x1j) Nx1
// subtract all rows of X from current row, col 1 only
r1 = (*CE.kernel)((I1[i] :- I1):/h[1])
dr1 = (*CE.dkernel)((I1[i] :- I1):/h[1]) / h[1]
// Leave one out
r1[i] =0
dr1[i]=0
// calculate kernel weights of other indeces, if any
if (args()==6) {
// r (X2i - X2j) Nxcols(X)-1
r2 = (*CE.kernel)((I2[i,.] :- I2) :/ h[(2..cols(h))] )
// Accumulate product across columns
for (k=1 ; k<=(cols(I2)) ; k++) {
r1 = r1:*r2[.,k]
dr1 = dr1:*r2[.,k]
}
}
// cexp[i,.] 1xcols(y)
real scalar dp , ydp
p[i] = quadsum(r1)
dp = quadcolsum(dr1)
yp[i] = cross(y,r1)
ydp = cross(y,dr1)
dcexpdI[i] = (ydp:*p[i] - yp[i]:*dp) / p[i]^2
}
}
// Add derivatives of all Lj wrt Ii for each row
for (i=1 ; i<=rows(y) ; i++) {
if (tx[i]) {
dKj = -(*CE.dkernel)((I1 :- I1[i]) /h[1]) /h[1]
if (args()==6) {
for (k=1 ; k<=(cols(I2)) ; k++) {
dKj = dKj :* (*CE.kernel)((I2[.,k] :- I2[i,k]) / h[k+1])
}
}
dLj = dKj :* (p*y[i] - yp) :/ p:^2
dLj[i]=0
dcexpdI[i] = dcexpdI[i] + quadsum(dLj)
}
}
return(dcexpdI)
}
/*******************************************************************************
Derivative of Conditional Expectation with respect to bandwidth
Vector Method
*******************************************************************************/
real matrix dcexpdh_vec(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
real vector h, dcexpdh
real scalar j
// Get bandwidth, row vector size K.
h = CE.optbwidth
// Initialize return matrix
dcexpdh = J(rows(X),cols(X),.)
// Loop through each column of h and X
for (j=1 ; j<=cols(h) ; j++) {
// Calculate derivative w.r.t. first bandwidth in vector
dcexpdh[.,j] = dcexpdh1_vec(y , X , CE , tx , h)
// Iterate throuch columns of X and h
// Move first column to last column, second to first, etc.
X = (X[.,(2..cols(X))],X[.,1])
h = (h[.,(2..cols(h))],h[.,1])
}
return(dcexpdh)
}
real vector dcexpdh1_vec(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx, real vector h) {
// Note: In this function X contains indexes
// I1 is the distance between two observations (the argument for the kernel function)
real scalar i , k
real vector r1, dr1, I1
real matrix dh1cexp , r
// Initialize cexp accumulation vector:
// value=0 for obs that will be trimmed.
dh1cexp = J(rows(X),cols(y),0)
// Calculate conditional expectation row by row.
for (i=1 ; i<=rows(X) ; i++) {
// Skip rows that will be trimmed.
if (tx[i]) {
// r1 (x1i - x1j) Nx1
// subtract all rows of X from current row, col 1 only
I1 = (X[i,1] :- X[.,1]):/h[1]
r1 = (*CE.kernel)(I1)
dr1 = (*CE.dkernel)(I1) :* -I1 :/ (h[1])
// note: I1 includes division by h, so no squared term needed in derivative
// Leave one out
r1[i] =0
dr1[i]=0
// calculate kernel weights of other X variables
if (cols(X)>1) {
// r (X2i - X2j) Nxcols(X)-1
r = (*CE.kernel)( ( X[i,(2..cols(X))] :- X[.,(2..cols(X))] ) :/ h[(2..cols(X))] )
// Accumulate product across columns
for (k=1 ; k<=(cols(X)-1) ; k++) {
r1 = r1:*r[.,k]
dr1 = dr1:*r[.,k]
}
}
// cexp[i,.] 1xcols(y)
real scalar p, dp
real rowvector ydp, yp
p = quadsum(r1)
dp = quadsum(dr1)
yp = cross(r1,y)
ydp = cross(dr1,y)
dh1cexp[i,.] = (ydp*p - yp*dp) / p^2
}
}
return(dh1cexp)
}
/*******************************************************************************
Kernel Functions
Note: Step-wise twicing procedure is as follows:
1. estimate E[y|x]
2. calculate r = y - E[y|x]
3. estimate E[r|x]
4. revise E*[y|x] = E[y|x]+E[r|x]
*******************************************************************************/
real kernel_twicing(real X) {
return(2*normalden(X) - normalden(X,0,sqrt(2)))
}
real dkernel_twicing(real X) {
return(-2*X :* normalden(X) + X/2 :* normalden(X,0,sqrt(2)))
}
real kernel_gaussian(real X) {
return(normalden(X))
}
real dkernel_gaussian(real X) {
return(-X:*normalden(X))
}
real kernel_epanechnikov(real X) {
real w
w = abs(X) :< 1
return(w :* 0.75 :* (1 :- X:^2))
}
real kernel_biweight(real X) {
real w
w = abs(X) :< 1
return(w :* 0.9375 :* (1 :- X:^2):^2)
}
real kernel_cosine(real X) {
real scalar twopi
real w
w = abs(X) :< 0.5
return(w :* cos(2*pi():*X))
}
real kernel_rectangular(real X) {
real w
w = abs(X) :< 1
return(w :* 0.5)
}
real kernel_triangular(real X) {
real absx, w
absx = abs(X)
w = absx :< 1
return(w :* (1-absx))
}
/*******************************************************************************
Bandwidth Functions
Constants from Hansen, Bruce E., Lecture Notes on Nonparametrics No. 1
To include higher order kernels in future update, kernels should be
specified by fxn and order.
Let k=num vars, and v=order
Then C should be specified in a v by k matrix.
and exponent should be n^(-1/(2*v+k))
*******************************************************************************/
real vector bwidth_pilot(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
// Declarations
real scalar n, k
real vector c, h, r
pointer scalar ptwicing, pgaussian, pepanech, pbiweight
n=rows(X)
k=cols(X)
// Get address of kernel fxns
ptwicing=&kernel_twicing()
// Choose c depending on kernel fxn and order
if ((CE.kernel)==ptwicing) {
// r ave values from KV 2010, p.158 (D1)
r = (.12916667 , .09166667)
h = sqrt(diagonal(variance(X))) :* n^(-r[k])
// c from Bruce for 2nd order gaussian
// c = (1.08,1.12)
// h = c[k] :* sqrt(diagonal(variance(X))) :* n^(-1/(8+k))
}
else {
pgaussian =&kernel_gaussian()
pepanech =&kernel_epanechnikov()
pbiweight =&kernel_biweight()
// ptriweight=&kernel_triweight()
if (CE.kernel==pgaussian ) c=(1.06,1.00)
else if (CE.kernel==pepanech ) c=(2.34,2.20)
else if (CE.kernel==pbiweight ) c=(2.78,2.61)
// else if (CE.kernel==ptriweight) c=(3.15,2.96)
else c=(1,1)
h = c[k] :* sqrt(diagonal(variance(X))) :* n^(-1/(4+k))
}
return(h')
}
/*******************************************************************************
SLS Bandwidth
Find optimal bandwidth (in terms of minimizing sum of squares) for each
set of potential parameters.
Estimate the log of h, rather than h.
This ensures that h can never be zero or negative.
*******************************************************************************/
void bwidth_sls(todo, real vector lnh, real vector y, real matrix X,
struct cexp_parameters scalar CE, real vector tx,
real vector minh, ssq, g, H) {
real vector h, est, sqdiff, trim, dEdh
// Set minimum bandwidth
h = exp(lnh) + minh
// Estimate Conditional Expectation and squared difference
est = (*CE.method)(y,X,CE,tx)
sqdiff = (y-est):^2
// Check for missing values
trim = missing(sqdiff)
if (trim>0) {
printf("%f observations trimmed due to zero denominator in cond. exp. \n", trim)
}
_editmissing(sqdiff, 0)
// Return objective fxn
ssq = tx :* sqdiff
if (todo==1) {
dEdh = dh1cexp(y,X,CE,tx,h)
g = -2 :* tx :* (y-est) :* dEdh
}
}
real vector bwidth_min_sls(real vector y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
real vector h, minh, lnh, delta
transmorphic scalar S
h=(bwidth_pilot(y,X,CE,tx))
minh = h :* 0.1
// db "Estimating Optimal Bandwidth"
// db "Starting value = "
// db h
// db "Minimum Bandwidth (1/10 pilot bandwidth)"
// db minh
lnh=ln(h-minh)
delta=J(1,cols(X),0.1)
S=optimize_init()
optimize_init_iterid(S, "Bandwidth")
optimize_init_valueid(S, "SSq(h)")
optimize_init_evaluatortype(S, "gf0")
optimize_init_which(S , "min")
optimize_init_evaluator(S, &bwidth_sls())
optimize_init_params(S, lnh)
optimize_init_tracelevel(S, "none")
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_argument(S, 3, CE)
optimize_init_argument(S, 4, tx)
optimize_init_argument(S, 5, minh)
// optimize_init_technique(S,"nr")
optimize_init_technique(S,"nm")
optimize_init_nmsimplexdeltas(S,delta)
// optimize_init_technique(S,"bhhh")
optimize_init_conv_ignorenrtol(S, "on")
lnh=optimize(S)
h = exp(lnh) + minh
// ptol default 1e-6
if (mreldif(h, minh) < optimize_init_conv_ptol(S)) {
printf("Bandwidth selection stopped at minimum (1/10 of pilot bw): \n")
h
}
// db "Optimal Bandwidth:"
// db h
return(h)
}
/*******************************************************************************
bandwidth Selection for derivative estimate.
Minimize score vector of sls objective fxn
*******************************************************************************/
void d1bwidth_sls(todo, real vector lnh, real vector y, real matrix X,
struct cexp_parameters scalar CE, real vector tx,
real vector minh, real vector r, ssq, g, H) {
real vector h, d1est, sqscore, trim
// Set minimum bandwidth
h = exp(lnh) + minh
CE.optbwidth = h
// Estimate Conditional Expectation and squared difference
d1est = d1cexp_vec(y,X,CE,tx)
sqscore = (r :* d1est):^2
// Check for missing values
trim = missing(sqscore)
if (trim>0) {
printf("%f observations trimmed due to zero denominator in cond. exp. \n", trim)
}
_editmissing(sqscore, 0)
// Return objective fxn
ssq = tx :* sqscore
}
real vector d1bwidth_min_sls(real vector y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
real vector h, minh, lnh, delta, Ey, r
transmorphic scalar S
// Expected value of Y and residual
Ey = cexp(y, X, CE, tx)
r = y-Ey
// Minimum bandwidth
h=(bwidth_pilot(y,X,CE,tx))
minh = h :* 0.1
// db "Estimating Optimal Bandwidth"
// db "Starting value = "
// db h
// db "Minimum Bandwidth (1/10 pilot bandwidth)"
// db minh
lnh=ln(h-minh)
delta=0.1
S=optimize_init()
optimize_init_iterid(S, "Derivative Bandwidth")
optimize_init_valueid(S, "SSq(h)")
optimize_init_evaluatortype(S, "gf0")
optimize_init_which(S , "min")
optimize_init_evaluator(S, &d1bwidth_sls())
optimize_init_params(S, lnh)
optimize_init_tracelevel(S, "step")
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_argument(S, 3, CE)
optimize_init_argument(S, 4, tx)
optimize_init_argument(S, 5, minh)
optimize_init_argument(S, 6, r)
// optimize_init_technique(S,"nr")
optimize_init_technique(S,"nm")
optimize_init_nmsimplexdeltas(S,delta)
optimize_init_conv_ignorenrtol(S, "on")
lnh=optimize(S)
h = exp(lnh) + minh
// ptol default 1e-6
if (mreldif(h, minh) < optimize_init_conv_ptol(S)) {
printf("Derivative Bandwidth selection stopped at minimum (1/10 of pilot bw): \n")
h
}
// db "Optimal Bandwidth:"
// db h
"Optimal Derivative Bandwidth:"
h
return(h)
}
/*******************************************************************************
Passthrough bandwidth.
Function to return current bandwidth parameter.
This can be used if the user wants to choose the bandwidth parameter
or if it is set by an outside optimization routine.
*******************************************************************************/
real vector bwidth_passthru(real matrix y , real matrix X ,
struct cexp_parameters scalar CE ,
real vector tx) {
// get current bandwidth vector
real vector h
h = CE.optbwidth
// ensure that a valid value for h exists.
if (length(h)==0) {
_error(3498, "optbwidth must be set before calling bwidth_passthru")
}
return(h)
}
/*******************************************************************************
Smooth Trimming
Return trimming vector
Multiply returned result by original data (element-wise) to get trimmed data
For large negative values in x, trimming weight is too small to be
represented as a double. More precisely, exp(-a*x) is larger than
maxdouble(), which results in a missing value.
In this case, trim to zero.
Note: Do not use smallestdouble(), which is smallest full precision
number. There are smaller, denormalized numbers that could exist for
trimming weights. I should not assign smallestdouble to missing values,
because it is not actually the smallest trimming weight available.
*******************************************************************************/
real vector strim(real vector x) {
real scalar n, a
real vector t
n = rows(x)
a = ln(n):^2
t = (1:+exp(-a:*x)):^(-1):*sign(x)
_editmissing(t,0)
return(t)
}
/*******************************************************************************
Indicator Trimming
Trim by centile. If any variable lies outside the specified centile range
for a given observation, that observation is omitted from the value function
calculation.
Returns a binary trimming vector.
I have to update this fxn for cases of repeated X variables. Currently,
this function identifies outliers by position, but if the same value lies
above and below the cut off position, they should be treated the same. Now,
one would be trimmed and one not trimmed.
Currently, indicator trimming is based on Stata's _pcentile function.
*******************************************************************************/
real vector itrim(real matrix X, real scalar lower_centile, real scalar upper_centile ) {
real scalar n, k, bn, tn, mn, j
real vector t, tx, invp
n=rows(X)
k=cols(X)
// number of obs to trim from bottom and top of distribution
bn = floor(lower_centile*n/100)
tn = floor((100-upper_centile)*n/100)
// number of untrimmed obs in the middle
mn = n-bn-tn
// ordered trimming vector
t = (J(bn,1,0) \ J(mn,1,1) \ J(tn,1,0))
// initialize accumulation trimming vector
tx = J(n,1,1)
// for each X
for (j=1 ; j<=k ; j++) {
// get inverse permutation vector...the vector to go from the ordered
// trimming vector to the current order of X[.,j]
invp = order(order(X,j),1)
// note: I could use the invorder function above, I think.
// Update accumulation vector with current vector, reordered to match X[j]
tx = tx :* t[invp,.]
}
return(tx)
}
/*******************************************************************************
Wrapper for calling itrim with variable names as strings.
_st_itrim returns a trimming vector in Mata
st_itrim generates a new variable in Stata
*******************************************************************************/
real vector _st_itrim(string scalar varlist, string scalar touse,
real scalar lower_centile, real scalar upper_centile) {
real vector tx
real matrix X
// Get data from Stata
st_view(X , . , varlist , touse)
tx = itrim(X,lower_centile,upper_centile)
st_numscalar("r(N)", rows(tx))
st_numscalar("r(sum)", sum(tx))
st_numscalar("r(mean)", mean(tx)[1] )
return(tx)
}
void st_itrim(string scalar varlist, string scalar touse,
real scalar lower_centile, real scalar upper_centile,
string scalar generate) {
real vector tx, idx
tx = _st_itrim(varlist, touse, lower_centile, upper_centile)
// Save trimming vector as new var
idx = st_addvar("byte", generate)
st_store(., idx , touse, tx)
}
/*******************************************************************************
Save function library
*******************************************************************************/
// mata mlib create lsls, dir(PERSONAL) replace
mata mlib create lsls, replace
mata mlib add lsls *()
// If no errors, 31 new functions were added.
mata mlib index
mata set matastrict off
end