Advice on replicating lpoly default bandwidth

[neal maroney]

Having difficulty replicating default bandwidth for lpoly with STATA R manual pg 1369.
It is an a function that is not that intuitive and not well described so i wanted to give it a shot.
can't quite get the right values.

Here are the steps (done in the attached do file) as I see them

1. polynomial degree =1
2. set constant term C01(Kernel=Guassian)= 0.76865817
3. run p+3 order global polynomial regression
4. calculate derivative
   der= (p + 1)! b[[p + 1]] + (p + 2)! b[[p +2]] x+ (p + 3)! b[[p + 3]] x^2/2
   on eq 5.4 Fan Gijbels (1995) verified with R locpol package. (STATA manual
   should document this step) Then sum the der^2 for denominator
5. get RSS from regression and std RSS of Schmidt[1973] for numerator candidate
6. ROT =c01K*(num/denom))^(1/(2*p + 3))
7. try on another dataset

I get in the neighborhood but cannot seem to get it right. Maybe the trimming of data issue - I don't know.
any advise?
thanks in advance

[Dirk Enzmann]

I can't answer your question, but: Please don't attach files! You can copy and paste syntax commands into your post and wrap [CODE] tags around the text of your syntax commands, see editor symbol #. See also the StataList FAQ, section 12.3.

[neal maroney]

K here is the code
thanks

* Neal Maroney 2/9/23 */
/* do file illustrating the replicating  of stata LPOLY default bandwidth HOPT */
/* try two types of data */
sca datafile="1";

if datafile =="1" {
    sca datais="  the datatset is motorcycle"
    use http://www.stata-press.com/data/r13/motorcycle, clear
    g y=accel /* ys */
    g x=time /*xs */
}
else {
    sca datais="   the datatset is fish"
    use https://www.stata-press.com/data/r17/fishdata, clear
    g y=freq /* ys */
    g x=length /*xs */
}


/*
Lpoly HOPT replication:
Refs
1. Fan and Gijbels "adaptive order Polynomial Fitting: Bandwidth Robustification
   and Bias Reduction", J. Computational Graphical Stats 1995. Eq 5.4 more precisely stated than eq
   in stata manual.
2. STAT lpoly R Manual p1369.  A bit vague on the details
3. Fan and Gijbels 1996
4. R Locpol package to verify derivatives. The R-locpol pluginBW[] is not the same as STATA even
   though the same resources are cited.  Different process entirely    


sca p=1
sca cvp01= 0.76865817 /*constant term for v=0,p=1 gaussian kernel*/

qui lpoly y x, kernel(gaussian) degree(1) nograph
sca HOPTLPOLY=r(bwidth)

/*global polynomial of order p+3 */
egen mx=mean(x)
sca dm=0 /* can demean dm=1 but didn't work*/

g o=1
g x1=(x-mx*dm)
g x2=(x-mx*dm)^2
g x3=(x-mx*dm)^3
g x4=(x-mx*dm)^4
qui reg y x1 x2 x3 x4
mat b= e(b)

predict e, residuals
mkmat e, matrix(es)
mat s2me=es'*es/(_N-5) /* regression variance */
sca  s2e=s2me[1,1]

/* standardized residual sum of squares Schmidt[1973]*/
mkmat o x1 x2 x3 x4, matrix(xs)
mat D=diag(vecdiag(I(_N)-xs*inv(xs'*xs)*xs')')
mat s2m=es'*inv(D)*es/(_N-5)
sca s2=s2m[1,1]

/*factorials needed*/
sca f2= 2*1  /* (p + 1)! */
sca f3= f2*3 /* (p + 2)! */
sca f4= f3*4 /* (p + 3)! */

/*m(p+1) derivative*/
g der =f2*  b[1,p+1] + f3 * b[1,p+2]*x + f4* b[1,p+3]*x^2/2
g der2=der^2 /* calls for summing the squared derivative*/
qui sum der2
sca sder2 =r(sum)

sca HOPT_SRRS=cvp01*(s2/sder2)^(1/(2*p + 3))
sca HOPT_REGVAR=cvp01*(s2e/sder2)^(1/(2*p + 3))

di datais
di  s2e   /* Regression variance- numerator */
di  s2    /* Std RSS/(n-k) -numerator*/
di  sder2 /* sum m(p+1) derivative squared- denominator */
dis HOPTLPOLY /*HOPT=Stata Lpoly */
dis HOPT_SRRS   /*HOPT=based on standardized rss. */
dis HOPT_REGVAR /*HOPT=based on regression variance */

/*
fish data results
HOPTLPOLY= 1.8414141
HOPT_SRRS  = 1.772881
HOPT_REGVAR =1.717532
*/

/*
motorcycle data results
HOPTLPOLY= 1.5445205
HOPT_SRRS  = 1.5537
HOPT_REGVAR =1.5434
*/

[neal maroney]

sorry, post in duplicate

[neal maroney]

How is the Standardized Residual Sum of Squares calculated ?
Stata R manual p1369.

Only reference for the equation I found is in
Hocking R.R., 1976, "The Analysis and Selection of Variables in Linear Regression", Biometrics, V32(1) p1-49.
https://www.jstor.org/stable/2529336
p16 cites Schmidt[1973] for standardized residual sum of squares:

srss = e'inv(D)*e, where D=DIAG(I - X*inv(X'X)*X')

I illustrate below for p=1.

use http://www.stata-press.com/data/r13/motorcycle, clear
g y=accel
g x=time

g o=1
g x1=x
g x2=x^2
g x3=x^3
g x4=x^4
qui reg y x1 x2 x3 x4
predict e, residuals
mkmat e, matrix(es)
mkmat o x1 x2 x3 x4, matrix(xs)

/*standardized residual sum of squares */

mat D=diag(vecdiag(I(_N)-xs*inv(xs'*xs)*xs')')
mat SRSS=es'*inv(D)*es
matlist SRSS /* = 213408.3 */
di e(rss) / *=206378.7 */

Results in slight difference in weighting compared to ordinary version.
Anyone know if this is the Stata calculation?