*! Source of lmoremata14.mlib *! {smcl} *! {marker mm_loclin}{bf:mm_loclin.mata}{asis} *! version 1.0.1 16jun2025 Ben Jann version 14.2 mata: real vector mm_loclin(real colvector Y, real colvector X, real colvector w, real vector at, real vector bw, | string scalar kernel) { real scalar i, varbw string scalar kern real vector fit pointer scalar k kern = _mm_unabkern(kernel) // default is "epan2" varbw = length(bw)!=1 fit = J(rows(at), cols(at), .) if (kern=="gaussian") { for (i=length(fit); i; i--) { fit[i] = _mm_loclin_gaussian(Y, X, w, at[i], varbw ? bw[i] : bw) } return(fit) } k = _mm_findkern(kern) for (i=length(fit); i; i--) { fit[i] = _mm_loclin(Y, X, w, at[i], varbw ? bw[i] : bw, k) } return(fit) } real scalar _mm_loclin_gaussian(real colvector Y, real colvector X, real colvector w, real scalar at, real scalar bw) { real colvector Z Z = X :- at return(mm_lsfit(Y, Z, w :* mm_kern_gaussian(Z/bw))[2]) } real scalar _mm_loclin(real colvector Y, real colvector X, real colvector w, real scalar at, real scalar bw, pointer scalar kern) { real scalar n real colvector W, Z, p Z = X :- at W = (*kern)(Z/bw) p = selectindex(W) n = length(p) if (!n) return(.) // no data within window if (n==rows(Y)) return(mm_lsfit(Y, Z, W:*w)[2]) return(mm_lsfit(Y[p], Z[p], (W:*w)[p])[2]) } real scalar mm_loclin_bw(real colvector Y, real colvector X, | real colvector w, string scalar kernel, real scalar fw) { real scalar h, c string scalar kern pointer scalar k if (args()<3) w = 1 if (args()<5) fw = 0 // compute canonical bandwidth of kernel; using mm_kdel0(kernel) would be // more precise and more efficient, but the goal is to reproduce lpoly's // bandwidth as closely as possible; thanks to Jeff Pitblado from Stata // Corp for pointing me to the existence of (undocumented) function // _lpoly_CnupK() kern = _mm_unabkern(kernel) // default is "epan2" if (kern=="epanechnikov") {; h = sqrt(5); k = &_lpoly_epanechnikov(); } else if (kern=="epan2") {; h = 1.25; k = &_lpoly_epan2(); } else if (kern=="biweight") {; h = 1.25; k = &_lpoly_biweight(); } else if (kern=="cosine") {; h = 1.25; k = &_lpoly_cosine(); } else if (kern=="gaussian") {; h = 5; k = &_lpoly_gaussian(); } else if (kern=="parzen") {; h = 1.25; k = &_lpoly_parzen(); } else if (kern=="rectangle") {; h = 1.25; k = &_lpoly_rectangle(); } else if (kern=="triangle") {; h = 1.25; k = &_lpoly_triangle(); } else _error(3300) // (cannot be reached) c = _lpoly_CnupK(-h, h, 2*h/100, 1, 0, k) // compute ROT bandwiath return(_mm_loclin_bw(Y, X, w, fw, 0) * c) } real scalar mm_loclin_bw2(real colvector Y, real colvector X, | real colvector w, string scalar kernel, real scalar fw) { if (args()<3) w = 1 if (args()<5) fw = 0 return(_mm_loclin_bw(Y, X, w, fw, 1) * mm_kdel0(kernel)) } real scalar _mm_loclin_bw(real colvector Y, real colvector X, real colvector w, real scalar fw, real scalar alt) { real scalar ll, ul, s2, M2, n real rowvector minmax real colvector b transmorphic S minmax = minmax(X) ll = minmax[1] + .05*(minmax[2] - minmax[1]) ul = minmax[2] - .05*(minmax[2] - minmax[1]) S = mm_ls(Y, (X,X:^2,X:^3,X:^4), w) b = mm_ls_b(S) s2 = mm_ls_rss(S)/mm_ls_N(S) n = fw ? (rows(w)==1 ? w*rows(X) : sum(w)) : rows(X) // Stata's lpoly takes the integral of the 2nd derivative of the global // polynomial fit without considering how X is distributed (i.e. assuming // an even distribution); however, according to the definition of the ROT // bandwidth by Fan and Gijbels (1996, page 111) the integral should be // taken across the distribution of X; this is what is done if alt!=0; else // the same approach as in Stata's lpoly is used if (alt) M2 = mean(_mm_loclin_bw_d2(X, b) :* (X:>=ll :& X:<=ul), w) else M2 = mm_integrate_sr(&_mm_loclin_bw_d2(), ll, ul, 1000, 0, b) return((s2 * (ul-ll) / (n * M2))^0.2) } real colvector _mm_loclin_bw_d2(real colvector X, real vector b) return((2*b[2] :+ 6*b[3]*X :+ 12*b[4]*X:^2):^2) end *! {smcl} *! {marker _mm_pieces14}{bf:_mm_pieces14.mata}{asis} *! version 1.0.0, Ben Jann, 01jun2015 version 14.0 mata: string rowvector _mm_pieces14( string scalar s, real scalar w, real scalar nobr) { real scalar i, j, k, n, l, b string scalar c string rowvector res l = udstrlen(s) if (w>=l | (l=ustrlen(s))<2) return(ustrtrim(s)) res = J(1, _mm_npieces14(s, w, nobr), "") j = k = n = 0 b = 1 for (i=1; i<=l; i++) { c = usubstr(s, i, 1) if (j<1) { // skip to first nonblank character if (ustrtrim(c)=="") { b++ continue } } j = j + udstrlen(c) if (i==l) { if (w>1 & !nobr) { // add extra row if last char is ud if (j>w) { res[++n] = ustrtrim(usubstr(s, b, i-b)) b = i } } res[++n] = ustrtrim(usubstr(s, b, .)) } else { if (ustrtrim(c)=="") k = i if (j>=w) { if (k<1) { if (nobr) continue if (i>b & j>w) k = i-1 else k = i } else { if (ustrtrim(usubstr(s, i+1, 1))=="") k = i } res[++n] = ustrtrim(usubstr(s, b, k-b+1)) j = udstrlen(usubstr(s, k+1, i-k)) b = k + 1 k = 0 } } } return(res) } real scalar _mm_npieces14( string scalar s, real scalar w, real scalar nobr) { real scalar i, j, k, n, l, b string scalar c l = udstrlen(s) if (w>=l | (l=ustrlen(s))<2) return(1) j = k = n = 0 b = 1 for (i=1; i<=l; i++) { c = usubstr(s, i, 1) if (j<1) { // skip to first nonblank character if (ustrtrim(c)=="") { b++ continue } } j = j + udstrlen(c) if (i==l) { if (w>1 & !nobr) { // add extra row if last char is ud if (j>w) n++ } n++ } else { if (ustrtrim(c)=="") k = i if (j>=w) { if (k<1) { if (nobr) continue if (i>b & j>w) k = i-1 else k = i } else { if (ustrtrim(usubstr(s, i+1, 1))=="") k = i } n++ j = udstrlen(usubstr(s, k+1, i-k)) b = k + 1 k = 0 } } } if (n==0) n++ return(n) } end