Return to Question

After inspecting the respective source codes I think the problem with Stata's implementation lies in fact that the subsample size used in each local regression depends on the ordinal position of the X values. That is, for extreme and near-extreme values of X (in other words, for values close to the tails), Stata uses a smaller subsample than for more central X's. Intuitively, the problem can be illustrated using a simple example with 100 data points where the bandwidth parameter is chosen to be 0.4 so that each subsample is of the size 0.4*100=40. In R, the size of the subsamples used for estimating $Y_1$, $Y_{10}$, $Y_{20}$, $Y_{30}$ and $Y_{40}$ would look like this:

On the other hand, Stata seems to reduce the size of the subsamples for all X's from (roughly) 1 to 20 and from 80 to 100:

When I rewrote the ksm program such that the subsample size was held fixed for all X's, I got the same results as in R (or Python). Below are the pseudo-codes of both implementations:

After inspecting the respective source codes I think the problem with Stata's implementation lies in fact that the subset size used in each local regression depends on the ordinal position of the X values. That is, for extreme and near-extreme values of X (in other words, for values close to the tails), Stata uses a smaller subset than for more central X's. Intuitively, the problem can be illustrated using a simple example with 100 data points where the bandwidth parameter is chosen to be 0.4 so that each subset is of the size 0.4*100=40. In R, the size of the subsets used for estimating $Y_1$, $Y_{10}$, $Y_{20}$, $Y_{30}$ and $Y_{40}$ would look like this:

On the other hand, Stata seems to reduce the size of the subsets for all X's from (roughly) 1 to 20 and from 80 to 100:

When I rewrote the ksm program such that the subsets size was held fixed for all X's, I got the same results as in R (or Python). Below are the pseudo-codes of both implementations:

Edit I: I wonder if this is because historically Stata's lowess command was based on the ksm program, which (when used with its default settings) was not meant to estimate the LOWESS regression. It appears that since Stata 8, the decision was made to abandon ksm and instead implement lowess where the default was chosen to be the LOWESS regression (whereas ksm's default options were implemented as lowess's non-defaults). The treatment of the tails, however, remained the same as in ksm.

Edit II: Note that all of these comparisons relied on matching all other parameters of the LOWESS regressions. As Stata doesn't allow for multiple iterations and also doesn't implement the interpolation to cut down on the number of local regressions required (as do R and the original Fortran code), I made sure R's command used both iter=0 and delta=0 parameters. Specifically, I was comparing lowess(X, Y, f=0.4, delta=0, iter=0) (R) with lowess Y X, bwidth(0.4) (Stata) and statsmodels.api.nonparametric.lowess(Y, X, frac=0.4, it=0, delta=0) (Python). Both X and Y were well-behaved.

EDIT I: I wonder if this is because historically Stata's lowess command was based on the ksm program, which (when used with its default settings) was not meant to estimate the LOWESS regression. It appears that since Stata 8, the decision was made to abandon ksm and instead implement lowess where the default was chosen to be the LOWESS regression (whereas ksm's default options were implemented as lowess's non-defaults). The treatment of the tails, however, remained the same as in ksm.

EDIT II: Note that all of these comparisons relied on matching all other parameters of the LOWESS regressions. As Stata doesn't allow for multiple iterations and also doesn't implement the interpolation to cut down on the number of local regressions required (as do R and the original Fortran code), I made sure R's command used both iter=0 and delta=0 parameters. Specifically, I was comparing lowess(X, Y, f=0.4, delta=0, iter=0) (R) with lowess Y X, bwidth(0.4) (Stata) and statsmodels.api.nonparametric.lowess(Y, X, frac=0.4, it=0, delta=0) (Python). Both X and Y were well-behaved.

UPDATE: OK, it seems I wasn't the first to notice this. I just found an old user-written adjksm command, which "is identical to ksm except that the bandwidth of the smoother is constant along the x-axis." See here. I also confirmed that adjksm Y X, bwidth(0.4) lowess generates the same output as R's lowess() using iter=0 and delta=0, up to the 4th decimal point.

I dove into the source code of the R's lowess() function (which seems to be based on Cleveland's original Fortran code found here: netlib (dot) org (slash) go (slash) lowess (dot) f) and the legacy Stata code for the ksm program found in the Stata 7 ado update file (found here: stata (dot) com (slash) support (slash) updates (slash) stata7 (slash) ado). Note that post-v7 Statas implemented the lowess command using the _LOWESS C routine that is not exposed to the user. I verified that running Stata 7's ksm command using the optional lowess argument in Stata 14 generates the same ("incorrect") output as running lowess directly. That is, ksm Y X, lowess is the same as lowess Y X.

I dove into the source code of the R's lowess() function (which seems to be based on Cleveland's original Fortran code found here) and the legacy Stata code for the ksm program found in the Stata 7 ado update file (found here). Note that post-v7 Statas implemented the lowess command using the _LOWESS C routine that is not exposed to the user. I verified that running Stata 7's ksm command using the optional lowess argument in Stata 14 generates the same ("incorrect") output as running lowess directly. That is, ksm Y X, lowess is the same as lowess Y X.