******************************************************************************** * LPDENSITY STATA PACKAGE -- lpdensity mata functions * Authors: Matias D. Cattaneo, Michael Jansson, Xinwei Ma ******************************************************************************** *!version 2.3.1 2022-05-03 ** do lpdensity_fn.ado capture mata mata drop lpdensity_fn() mata real matrix lpdensity_fn(real colvector data, real colvector grid, real colvector bw, real scalar p, real scalar q, real scalar v, string scalar kernel, real colvector Cweights, real colvector Pweights, real scalar massPoints, real scalar scale, real scalar size, real scalar CIuniform, real scalar CIsimul){ ii = order(data, 1) data = data[ii, 1] Cweights = Cweights[ii, 1] Pweights = Pweights[ii, 1] n = length(data) ng = length(grid) dataUnique = lpdensity_unique(data) freqUnique = dataUnique[., 2] indexUnique = dataUnique[., 3] dataUnique = dataUnique[., 1] nUnique = length(dataUnique) if (massPoints) { Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique) } else{ Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights) } weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n Cweights = Cweights :/ sum(Cweights) :* n Pweights = Pweights :/ sum(Pweights) :* n weights_normalUnique = runningsum(weights_normal)[indexUnique] if (nUnique > 1) { weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) } CweightsUnique = runningsum(Cweights)[indexUnique] if (nUnique > 1) { CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) } PweightsUnique = runningsum(Pweights)[indexUnique] if (nUnique > 1) { PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) } Result = J(ng, 10, .) iff_p = J(n, ng, .) // store influence functions iff_q = iff_p // loop over evaluation points for (j=1; j<=ng; j++) { Result[j, 10] = j Result[j, 1] = grid[j] Result[j, 2] = bw[j] index_temp = (abs(data :- grid[j]) :<= bw[j]) Result[j, 3] = sum(index_temp) // LHS and RHS variables if (massPoints) { Y_temp = Fn[indexUnique] Xh_temp = (data[indexUnique] :- grid[j]) :/ bw[j] // weights if (kernel == "triangular") { Kh_temp = ((1 :- abs(Xh_temp)) :/ bw[j]) :* index_temp[indexUnique] } if (kernel == "uniform") { Kh_temp = (0.5 :/ bw[j]) :* index_temp[indexUnique] } if (kernel == "epanechnikov") { Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ bw[j]) :* index_temp[indexUnique] } Xh_p_temp = J(nUnique, p+1, .) Xh_p_Kh_Pweights_temp = Xh_p_temp for (jj=1; jj<=p+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique } XhKhXh_inv = invsym(cross(Xh_p_temp, Xh_p_Kh_Pweights_temp) :/ n) // point estimate Result[j, 4] = factorial(v) :* (XhKhXh_inv * cross(Xh_p_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n // standard error estimate F_Xh_p_Kh_temp = cross(Y_temp, Xh_p_Kh_Pweights_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1; jj++) { G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal } iff_p[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))' } else { Y_temp = Fn Xh_temp = (data :- grid[j]) :/ bw[j] // weights if (kernel == "triangular") { Kh_temp = ((1 :- abs(Xh_temp)) :/ bw[j]) :* index_temp } if (kernel == "uniform") { Kh_temp = (0.5 :/ bw[j]) :* index_temp } if (kernel == "epanechnikov") { Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ bw[j]) :* index_temp } Xh_p_temp = J(n, p+1, .) Xh_p_Kh_Pweights_temp = Xh_p_temp for (jj=1; jj<=p+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights } XhKhXh_inv = invsym(cross(Xh_p_temp, Xh_p_Kh_Pweights_temp) :/ n) // point estimate Result[j, 4] = factorial(v) :* (XhKhXh_inv * cross(Xh_p_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n // standard error estimate F_Xh_p_Kh_temp = cross(Y_temp, Xh_p_Kh_Pweights_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1; jj++) { G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal } iff_p[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))' } // robust CI if (q > p) { if (massPoints) { Xh_q_temp = J(nUnique, q+1, .) Xh_q_Kh_Pweights_temp = Xh_q_temp for (jj=1; jj<=q+1; jj++) { Xh_q_temp[., jj] = Xh_temp:^(jj-1) Xh_q_Kh_Pweights_temp[., jj] = Xh_q_temp[., jj] :* Kh_temp :* PweightsUnique } XhKhXh_inv = invsym(cross(Xh_q_temp, Xh_q_Kh_Pweights_temp) :/ n) // point estimate Result[j, 5] = factorial(v) :* (XhKhXh_inv * cross(Xh_q_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n // standard error estimate F_Xh_q_Kh_temp = cross(Y_temp, Xh_q_Kh_Pweights_temp) :/ n G = J(n, q - 0 + 1, .) for (jj=1; jj<=q-0+1; jj++) { G[., jj] = (lpdensity_rep((runningsum(Xh_q_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_q_Kh_temp[jj]) :* weights_normal } iff_q[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))' } else { Xh_q_temp = J(n, q+1, .) Xh_q_Kh_Pweights_temp = Xh_q_temp for (jj=1; jj<=q+1; jj++) { Xh_q_temp[., jj] = Xh_temp:^(jj-1) Xh_q_Kh_Pweights_temp[., jj] = Xh_q_temp[., jj] :* Kh_temp :* Pweights } XhKhXh_inv = invsym(cross(Xh_q_temp, Xh_q_Kh_Pweights_temp) :/ n) // point estimate Result[j, 5] = factorial(v) :* (XhKhXh_inv * cross(Xh_q_Kh_Pweights_temp, Y_temp))[v+1] :/ bw[j]:^v :/ n // standard error estimate F_Xh_q_Kh_temp = cross(Y_temp, Xh_q_Kh_Pweights_temp) :/ n G = J(n, q - 0 + 1, .) for (jj=1; jj<=q-0+1; jj++) { G[., jj] = ((runningsum(Xh_q_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_q_Kh_temp[jj]) :* weights_normal } iff_q[., j] = ((XhKhXh_inv * G')[v+1, .] :* factorial(v) :/ sqrt(n * bw[j]^(2*v)))' } } } CovMat_p = cross(iff_p, iff_p) :/ n CovMat_q = cross(iff_q, iff_q) :/ n Result[., 6] = sqrt(abs(diagonal(CovMat_p))) Result[., 7] = sqrt(abs(diagonal(CovMat_q))) Result[., 4..7] = Result[., 4..7] :* scale if (CIuniform == 0) { critVal = invnormal(1-size/2) } else { if (q > p) { CovMat_q = CovMat_q :* scale^2 for (j1=1; j1 <= ng; j1++) { for (j2=1; j2 <= ng; j2++) { CovMat_q[j1, j2] = CovMat_q[j1, j2] / Result[j1, 7] / Result[j2, 7] } } normRVs = cholesky(CovMat_q) * rnormal(ng, CIsimul, 0, 1) critVal = lpdensity_quantile(colmax(abs(normRVs))', 1-size) } else { CovMat_p = CovMat_p :* scale^2 for (j1=1; j1 <= ng; j1++) { for (j2=1; j2 <= ng; j2++) { CovMat_p[j1, j2] = CovMat_p[j1, j2] / Result[j1, 6] / Result[j2, 6] } } normRVs = cholesky(CovMat_p) * rnormal(ng, CIsimul, 0, 1) critVal = lpdensity_quantile(colmax(abs(normRVs))', 1-size) } } if (q > p) { Result[., 8] = Result[., 5] - critVal * Result[., 7] Result[., 9] = Result[., 5] + critVal * Result[., 7] } if (q == p) { Result[., 8] = Result[., 4] - critVal * Result[., 6] Result[., 9] = Result[., 4] + critVal * Result[., 6] } return(Result) } mata mosave lpdensity_fn(), replace end ******************************************************************************** * Extracting unique elements ******************************************************************************** capture mata mata drop lpdensity_unique() mata real matrix lpdensity_unique(real colvector x){ // Note: x should be in ascending order n = length(x) // x has one or no element if (n == 0) { out = J(0, 3, .) return(out) } if (n == 1) { out = (x, 1, 1) return(out) } // else numIndexTemp = selectindex(x[2..n] :!= x[1..(n-1)]) if (length(numIndexTemp) == 0) { numIndexLast = J(1, 1, n) numIndexFirst = J(1, 1, 1) } else { numIndexLast = (numIndexTemp \ n) numIndexFirst = (1 \ numIndexTemp:+1) } freq = numIndexLast - numIndexFirst :+ 1 unique = x[numIndexLast] out = (unique, freq, numIndexLast) return(out) } /* // testing mata: x = (1\ 1) mata: lpdensity_unique(x) mata: x = (1\ 2\ 2) mata: lpdensity_unique(x) mata: x = (1..100)' mata: lpdensity_unique(x) mata: x = (1\ 1\ 2\ 3\ 4\ 4\ 5\ 6\ 6\ 6\ 7) mata: lpdensity_unique(x) mata: x = (1\ 1\ 2\ 3\ 4\ 4\ 5\ 6\ 6\ 6\ 7\ 7) mata: lpdensity_unique(x) mata: x = (5) mata: lpdensity_unique(x) */ mata mosave lpdensity_unique(), replace end ******************************************************************************** * Select the minimum element ******************************************************************************** capture mata mata drop lpdensity_whichmin() mata real scalar lpdensity_whichmin(real colvector x){ if (length(x) == 0) { return(0) } else { return (order(x, 1)[1]) } } /* // testing mata: x = (1\ 1) mata: lpdensity_whichmin(x) mata: x = (1\ 2\ 2) mata: lpdensity_whichmin(x) mata: x = (3\ 2\ 2) mata: lpdensity_whichmin(x) mata: x = (3\ 2\ 3) mata: lpdensity_whichmin(x) mata: x = (3\ 2\ 1) mata: lpdensity_whichmin(x) */ mata mosave lpdensity_whichmin(), replace end ******************************************************************************** * Replicating each element in x y times (elementwise) ******************************************************************************** capture mata mata drop lpdensity_rep() mata real matrix lpdensity_rep(real colvector x, real colvector y){ // Warning: x and y should have the same length // y should contain strictly positive integers // y cannot be empty nout = sum(y) if (length(y) == 1) { out = J(nout, 1, x) } else { if (all(y :== 1)) { out = x } else { out = J(nout, 1, .) out[1..y[1], 1] = J(y[1], 1, x[1]) indextemp = y[1] for (i=2; i<=length(y); i++) { out[(indextemp+1)..(indextemp+y[i])] = J(y[i], 1, x[i]) indextemp = indextemp + y[i] } } } return(out) } /* // testing mata: lpdensity_rep(5, 1) mata: lpdensity_rep(5, 3) mata: lpdensity_rep((5\ 6), (1\ 1)) mata: lpdensity_rep((5\ 6), (1\ 2)) mata: lpdensity_rep((5\ 6), (2\ 1)) mata: lpdensity_rep((5\ 6), (2\ 3)) mata: lpdensity_rep((5\ 6\ 7\ 8), (1\ 1\ 1\ 1)) mata: lpdensity_rep((5\ 6\ 7\ 8), (1\ 2\ 2\ 3)) mata: lpdensity_rep((5\ 6\ 7\ 8), (3\ 1\ 2\ 1)) mata: lpdensity_rep((5\ 6\ 7\ 8), (3\ 2\ 2\ 3)) */ mata mosave lpdensity_rep(), replace end ******************************************************************************** * Normal DGPs ******************************************************************************** capture mata mata drop lpdensity_norm_dgps() mata real scalar lpdensity_normDgp(real scalar y, real scalar v, real scalar mu, real scalar sigma){ x = (y - mu) / sigma // normalized evaluation point if (v == 0 ) out = normal(x) // the cdf if (v > 0) { // compute hermite polynomials if (v == 1) H = 1 if (v == 2) H = x if (v > 2) { temp = J(v, 1, .) temp[1] = 1; temp[2] = x for (j=3; j<=v; j++) temp[j] = x * temp[j-1] - (j-2) * temp[j-2] H = temp[v] } // get the sign correct H = H * (-1)^(v + 1) // get the scaling correct out = H * normalden(x) / sigma^v } return(out) } mata mosave lpdensity_normDgp(), replace end ******************************************************************************** * empirical quantile ******************************************************************************** capture mata mata drop lpdensity_quantile() mata real scalar lpdensity_quantile(real colvector x, real scalar p){ x = sort(x, 1) n = length(x) q = ceil(p * n) if (q < 1) q = 1 if (q > n) q = n out = x[q] return(out) } mata mosave lpdensity_quantile(), replace end ******************************************************************************** * S matrix generation ******************************************************************************** capture mata mata drop lpdensity_Sgen() mata real matrix lpdensity_Sgen(real scalar p, string scalar kernel){ if (kernel == "uniform") { out = (1,0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909\ 0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0\ 0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769\ 0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0\ 0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667\ 0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0\ 0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647\ 0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0\ 0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684\ 0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0\ 0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476) } if (kernel == "triangular") { out = (1,0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831\ 0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0\ 0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887\ 0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0\ 0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318\ 0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0\ 0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977\ 0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0\ 0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991\ 0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0\ 0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433) } if (kernel == "epanechnikov") { out = (1,0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021\ 0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0\ 0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154\ 0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0\ 0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529\ 0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0\ 0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443\ 0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0\ 0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812\ 0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0\ 0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236) } return(out[1..(p+1), 1..(p+1)]) } mata mosave lpdensity_Sgen(), replace end ******************************************************************************** * T matrix generation ******************************************************************************** capture mata mata drop lpdensity_Tgen() mata real matrix lpdensity_Tgen(real scalar p, string scalar kernel){ if (kernel == "uniform") { out = (0.5,0,0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455\ 0,0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0\ 0.166666666666667,0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385\ 0,0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0\ 0.1,0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333\ 0,0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0\ 0.0714285714285714,0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824\ 0,0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0\ 0.0555555555555556,0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842\ 0,0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842,0\ 0.0454545454545455,0,0.0384615384615385,0,0.0333333333333333,0,0.0294117647058824,0,0.0263157894736842,0,0.0238095238095238) } if (kernel == "triangular") { out = (0.666666666666667,0,0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814\ 0,0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0\ 0.0666666666666667,0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681\ 0,0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0\ 0.019047619047619,0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244\ 0,0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0\ 0.00793650347516803,0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965\ 0,0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0\ 0.00404040423366886,0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208\ 0,0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208,0\ 0.00233100231773814,0,0.00146520146655681,0,0.00098039215666244,0,0.000687994496086965,0,0.00050125313283208,0,0.000376435159043855) } if (kernel == "epanechnikov") { out = (0.6,0,0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042\ 0,0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0\ 0.0857142857142857,0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683\ 0,0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0\ 0.0285714285714286,0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889\ 0,0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0\ 0.012987012987013,0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492\ 0,0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0\ 0.00699300699300699,0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932\ 0,0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932,0\ 0.0041958041958042,0,0.00271493212669683,0,0.00185758513931889,0,0.00132684652808492,0,0.0009807126511932,0,0.000745341614906832) } return(out[1..(p+1), 1..(p+1)]) } mata mosave lpdensity_Tgen(), replace end ******************************************************************************** * C matrix generation ******************************************************************************** capture mata mata drop lpdensity_Cgen() mata real matrix lpdensity_Cgen(real scalar k, real scalar p, string scalar kernel){ if (kernel == "uniform") { out = (1,0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0\ 0,0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647\ 0.333333333333333,0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0\ 0,0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684\ 0.2,0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0\ 0,0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476\ 0.142857142857143,0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0\ 0,0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652\ 0.111111111111111,0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0\ 0,0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0,0.04\ 0.0909090909090909,0,0.0769230769230769,0,0.0666666666666667,0,0.0588235294117647,0,0.0526315789473684,0,0.0476190476190476,0,0.0434782608695652,0,0.04,0) } if (kernel == "triangular") { out = (1,0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0\ 0,0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977\ 0.166666666666667,0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0\ 0,0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991\ 0.0666666666666667,0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0\ 0,0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433\ 0.0357142834836158,0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0\ 0,0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971\ 0.0222222223188546,0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0\ 0,0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0,0.00307692307692308\ 0.0151515151448831,0,0.0109890109896887,0,0.00833333333323318,0,0.00653594771243977,0,0.00526315789472991,0,0.00432900432900433,0,0.0036231884057971,0,0.00307692307692308,0) } if (kernel == "epanechnikov") { out = (1,0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0\ 0,0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443\ 0.2,0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0\ 0,0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812\ 0.0857142857142857,0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0\ 0,0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236\ 0.0476190476190476,0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0\ 0,0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783\ 0.0303030303030303,0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0\ 0,0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0,0.00444444444444445\ 0.020979020979021,0,0.0153846153846154,0,0.0117647058823529,0,0.00928792569659443,0,0.0075187969924812,0,0.0062111801242236,0,0.00521739130434783,0,0.00444444444444445,0) } return(out[1..(p+1), k+1]) } mata mosave lpdensity_Cgen(), replace end ******************************************************************************** * G matrix generation ******************************************************************************** capture mata mata drop lpdensity_Ggen() mata real matrix lpdensity_Ggen(real scalar p, string scalar kernel){ if (kernel == "uniform") { out = (-0.333333333333333,0.166666666666667,-0.133333333333333,0.1,-0.0857142857142857,0.0714285714285714,-0.0634920634920635,0.0555555555555556,-0.0505050505050505,0.0454545454545455,-0.041958041958042\ 0.166666666666667,0.0666666666666667,0.0555555555555556,0.0380952380952381,0.0333333333333333,0.0264550264550265,0.0238095238095238,0.0202020202020202,0.0185185185185185,0.0163170163170163,0.0151515151515152\ -0.133333333333333,0.0555555555555556,-0.0476190476190476,0.0333333333333333,-0.0296296296296296,0.0238095238095238,-0.0216450216450216,0.0185185185185185,-0.0170940170940171,0.0151515151515152,-0.0141414141414141\ 0.1,0.0380952380952381,0.0333333333333333,0.0222222222222222,0.02,0.0155844155844156,0.0142857142857143,0.011965811965812,0.0111111111111111,0.0096969696969697,0.00909090909090909\ -0.0857142857142857,0.0333333333333333,-0.0296296296296296,0.02,-0.0181818181818182,0.0142857142857143,-0.0131868131868132,0.0111111111111111,-0.0103703703703704,0.00909090909090909,-0.00855614973262032\ 0.0714285714285714,0.0264550264550265,0.0238095238095238,0.0155844155844156,0.0142857142857143,0.010989010989011,0.0102040816326531,0.00846560846560847,0.00793650793650794,0.00687547746371276,0.0064935064935065\ -0.0634920634920635,0.0238095238095238,-0.0216450216450217,0.0142857142857143,-0.0131868131868132,0.0102040816326531,-0.00952380952380953,0.00793650793650794,-0.00746965452847806,0.0064935064935065,-0.00615174299384826\ 0.0555555555555556,0.0202020202020202,0.0185185185185185,0.011965811965812,0.0111111111111111,0.00846560846560847,0.00793650793650794,0.0065359477124183,0.00617283950617284,0.00531632110579479,0.00505050505050505\ -0.0505050505050505,0.0185185185185185,-0.0170940170940171,0.0111111111111111,-0.0103703703703704,0.00793650793650794,-0.00746965452847806,0.00617283950617284,-0.00584795321637427,0.00505050505050505,-0.00481000481000481\ 0.0454545454545455,0.0163170163170163,0.0151515151515152,0.0096969696969697,0.00909090909090909,0.00687547746371276,0.00649350649350649,0.00531632110579479,0.00505050505050505,0.00432900432900433,0.00413223140495868\ -0.041958041958042,0.0151515151515152,-0.0141414141414141,0.00909090909090909,-0.00855614973262032,0.00649350649350649,-0.00615174299384826,0.00505050505050505,-0.00481000481000481,0.00413223140495868,-0.00395256916996048) } if (kernel == "triangular") { out = (-0.233332866779462,0.0833330508905063,-0.0531743640785073,0.0333333960560459,-0.0243382259715177,0.0178568590991684,-0.0140329271889015,0.0111109568092491,-0.00914256712672596,0.00757561833195314,-0.00643232469904509\ 0.0833333430031619,0.0206349536543496,0.0138889822607069,0.00747358076709329,0.00555558002800799,0.00376385059786869,0.00297624095083857,0.00224683629513331,0.00185187157177693,0.00148741419069681,0.00126263242255087\ -0.0531747120633688,0.0138889436189969,-0.00992067054772658,0.00555558570938115,-0.00427612357529663,0.00297621392905644,-0.00240733368919841,0.00185187029011056,-0.00155068139981584,0.00126263712216438,-0.00108438068750666\ 0.0333333291761154,0.00747354338303244,0.0055555555888295,0.00282828659007681,0.00222222610846186,0.00145548770075321,0.00119047711034882,0.000879119696771778,0.000740741774432789,0.000586005979164481,0.00050505152540861\ -0.0243386235494463,0.00555555541211971,-0.00427609423962765,0.0022222251541454,-0.00178710387992154,0.00119047600313231,-0.000989883016025128,0.000740740651962696,-0.00063180824060111,0.000505050436046404,-0.000439378283083419\ 0.0178571428772336,0.0037638288080394,0.00297619046584015,0.00145549115212934,0.00119047780581393,0.000758765013497074,0.000637755044256656,0.000461990375582904,0.000396825380503709,0.000309560314164898,0.000270562752861001\ -0.014033189060362,0.00297619047992443,-0.0024073149074495,0.00119047782568061,-0.000989884403822643,0.000637755064921552,-0.00054271705691498,0.000396825400688911,-0.000344133756729751,0.000270562772022687,-0.000238285106683998\ 0.01111111111111,0.00224682724469665,0.00185185185231585,0.000879122147763556,0.000740741756863993,0.000461990319945608,0.000396825372642262,0.000282842182321384,0.000246913581918653,0.000190248347383829,0.000168350168854074\ -0.00914270914122569,0.00185185185168499,-0.00155067155074166,0.000740741756146379,-0.000631809109714714,0.000396825371926143,-0.000344133734564222,0.000246913581232753,-0.00021720969177002,0.000168350168212254,-0.000149908647886651\ 0.00757575757574178,0.00148740148756518,0.00126256257569384,0.000586007961071109,0.000505051197433514,0.000309560273844802,0.000270562753630222,0.000190248349061672,0.00016835016904926,0.000128330167831491,0.000114784205611466\ -0.00643245643257991,0.00126256257573434,-0.00108431925428119,0.000505051197475794,-0.000439378889047697,0.000270562753671301,-0.000238285092406906,0.000168350169087952,-0.000149908648550576,0.000114784205647198,-0.000103205972768449) } if (kernel == "epanechnikov") { out = (-0.257142857142857,0.1,-0.0666666666666667,0.0428571428571429,-0.032034632034632,0.0238095238095238,-0.018981018981019,0.0151515151515152,-0.0125874125874126,0.0104895104895105,-0.00896750308515014\ 0.1,0.0285714285714286,0.02,0.0112554112554113,0.00857142857142857,0.00592740592740593,0.00476190476190476,0.00363636363636364,0.00303030303030303,0.00245166598107775,0.0020979020979021\ -0.0666666666666667,0.02,-0.0147186147186147,0.00857142857142857,-0.00672660672660673,0.00476190476190476,-0.0039027639027639,0.00303030303030303,-0.00256136020841903,0.0020979020979021,-0.00181428478642101\ 0.0428571428571429,0.0112554112554113,0.00857142857142857,0.00459540459540459,0.0036734693877551,0.00246420246420246,0.00204081632653061,0.00152710035062976,0.0012987012987013,0.00103586815661119,0.000899100899100899\ -0.032034632034632,0.00857142857142857,-0.00672660672660672,0.0036734693877551,-0.002997002997003,0.00204081632653061,-0.00171514759750054,0.0012987012987013,-0.00111669548202056,0.000899100899100899,-0.000787094374002821\ 0.0238095238095238,0.00592740592740593,0.00476190476190476,0.00246420246420246,0.00204081632653061,0.00133591898297781,0.00113378684807256,0.000833634889362443,0.000721500721500722,0.000568059391588803,0.0004995004995005\ -0.018981018981019,0.00476190476190476,-0.0039027639027639,0.00204081632653061,-0.00171514759750054,0.00113378684807256,-0.000973020787262273,0.000721500721500722,-0.000629917038585769,0.0004995004995005,-0.000442294982994558\ 0.0151515151515152,0.00363636363636364,0.00303030303030303,0.00152710035062976,0.0012987012987013,0.000833634889362443,0.000721500721500721,0.000522697117124362,0.000459136822773186,0.000357384796744064,0.000317863954227591\ -0.0125874125874126,0.00303030303030303,-0.00256136020841903,0.0012987012987013,-0.00111669548202056,0.000721500721500721,-0.000629917038585769,0.000459136822773186,-0.000406153724231527,0.000317863954227591,-0.000284353327831589\ 0.0104895104895105,0.00245166598107775,0.0020979020979021,0.00103586815661119,0.000899100899100899,0.000568059391588803,0.000499500499500499,0.000357384796744064,0.000317863954227591,0.000244972418885462,0.000220059660619101\ -0.00896750308515014,0.0020979020979021,-0.00181428478642101,0.000899100899100899,-0.000787094374002822,0.000499500499500499,-0.000442294982994558,0.000317863954227591,-0.000284353327831589,0.000220059660619101,-0.000198641937772373) } return(out[1..(p+1), 1..(p+1)]) } mata mosave lpdensity_Ggen(), replace end ******************************************************************************** * mse-rot bandwidth ******************************************************************************** capture mata mata drop lpdensity_bwROT() mata real matrix lpdensity_bwROT(real colvector data, real colvector grid, real scalar p, real scalar v, string scalar kernel, real colvector Cweights, real colvector Pweights, real scalar massPoints, real scalar stdVar, real scalar regularize, real scalar nLocalMin, real scalar nUniqueMin){ data = sort(data, 1) n = length(data) ng = length(grid) dataUnique = lpdensity_unique(data)[., 1] nUnique = length(dataUnique) if (stdVar) { center_temp = mean(data) scale_temp = sqrt(variance(data)) data = (data :- center_temp) :/ scale_temp dataUnique = (dataUnique :- center_temp) :/ scale_temp grid = (grid :- center_temp) :/ scale_temp } // estimate a normal reference model mean_hat = sum(data :* Cweights :* Pweights) / sum(Cweights :* Pweights) sd_hat = sqrt(sum(Cweights :* Pweights :* (data :- mean_hat):^2) / sum(Cweights :* Pweights)) // bias estimate, no rate added bias_dgp = J(ng, 2, .) for (j=1; j<=ng; j++) { bias_dgp[j, 1] = lpdensity_normDgp(grid[j], p+1, mean_hat, sd_hat) / factorial(p+1) * factorial(v) bias_dgp[j, 2] = lpdensity_normDgp(grid[j], p+2, mean_hat, sd_hat) / factorial(p+2) * factorial(v) + bias_dgp[j, 1] * lpdensity_normDgp(grid[j], 2, mean_hat, sd_hat) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) } S = lpdensity_Sgen( p, kernel=kernel) C1 = lpdensity_Cgen(p+1, p, kernel=kernel) C2 = lpdensity_Cgen(p+2, p, kernel=kernel) G = lpdensity_Ggen( p, kernel=kernel) S2 = lpdensity_Tgen( p, kernel=kernel) bias_dgp[., 1] = bias_dgp[., 1] :* (invsym(S) * C1)[v+1, .] bias_dgp[., 2] = bias_dgp[., 2] :* (invsym(S) * C2)[v+1, .] // variance estimate, sample size added sd_dgp = J(ng, 1, .) if (v > 0) { for (j=1; j<=ng; j++) { sd_dgp[j, 1] = factorial(v) * sqrt(lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / n) } sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * G * invsym(S))[v+1, v+1])) } else { for (j=1; j<=ng; j++) { sd_dgp[j, 1] = factorial(v) * sqrt( lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat) * (1 - lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat)) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / (0.5 * n^2)) } sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * S2 * invsym(S))[v+1, v+1])) } // bandwidth h = J(ng, 1, .) for (j=1; j<=ng; j++) { S = optimize_init() optimize_init_evaluator(S, &lpdensity_optFunc()) optimize_init_params(S, 1e-8) optimize_init_argument(S, 1, bias_dgp[j, .]) optimize_init_argument(S, 2, sd_dgp[j, .]) optimize_init_argument(S, 3, p) optimize_init_argument(S, 4, v) optimize_init_argument(S, 5, "mse-rot") optimize_init_which(S, "min") optimize_init_evaluatortype(S, "gf0") optimize_init_tracelevel(S, "none") optimize_init_conv_maxiter(S, 20) (void) optimize(S) h[j] = optimize_result_params(S) if (regularize == 1) { if (nLocalMin > 0) { h[j] = max((h[j], sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))])) } if (nUniqueMin > 0) { h[j] = max((h[j], sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))])) } } h[j] = min(( h[j], max( abs(dataUnique :- grid[j]) ) )) } if (stdVar) { h = h :* scale_temp } return(h) } mata mosave lpdensity_bwROT(), replace end ******************************************************************************** * imse-rot bandwidth ******************************************************************************** capture mata mata drop lpdensity_bwIROT() mata real scalar lpdensity_bwIROT(real colvector data, real colvector grid, real scalar p, real scalar v, string scalar kernel, real colvector Cweights, real colvector Pweights, real scalar massPoints, real scalar stdVar, real scalar regularize, real scalar nLocalMin, real scalar nUniqueMin){ data = sort(data, 1) n = length(data) ng = length(grid) dataUnique = lpdensity_unique(data)[., 1] nUnique = length(dataUnique) if (stdVar) { center_temp = mean(data) scale_temp = sqrt(variance(data)) data = (data :- center_temp) :/ scale_temp dataUnique = (dataUnique :- center_temp) :/ scale_temp grid = (grid :- center_temp) :/ scale_temp } // estimate a normal reference model mean_hat = sum(data :* Cweights :* Pweights) / sum(Cweights :* Pweights) sd_hat = sqrt(sum(Cweights :* Pweights :* (data :- mean_hat):^2) / sum(Cweights :* Pweights)) // bias estimate, no rate added bias_dgp = J(ng, 2, .) for (j=1; j<=ng; j++) { bias_dgp[j, 1] = lpdensity_normDgp(grid[j], p+1, mean_hat, sd_hat) / factorial(p+1) * factorial(v) bias_dgp[j, 2] = lpdensity_normDgp(grid[j], p+2, mean_hat, sd_hat) / factorial(p+2) * factorial(v) + bias_dgp[j, 1] * lpdensity_normDgp(grid[j], 2, mean_hat, sd_hat) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) } S = lpdensity_Sgen( p, kernel=kernel) C1 = lpdensity_Cgen(p+1, p, kernel=kernel) C2 = lpdensity_Cgen(p+2, p, kernel=kernel) G = lpdensity_Ggen( p, kernel=kernel) S2 = lpdensity_Tgen( p, kernel=kernel) bias_dgp[., 1] = bias_dgp[., 1] :* (invsym(S) * C1)[v+1, .] bias_dgp[., 2] = bias_dgp[., 2] :* (invsym(S) * C2)[v+1, .] // variance estimate, sample size added sd_dgp = J(ng, 1, .) if (v > 0) { for (j=1; j<=ng; j++) { sd_dgp[j, 1] = factorial(v) * sqrt(lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / n) } sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * G * invsym(S))[v+1, v+1])) } else { for (j=1; j<=ng; j++) { sd_dgp[j, 1] = factorial(v) * sqrt( lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat) * (1 - lpdensity_normDgp(grid[j], 0, mean_hat, sd_hat)) / lpdensity_normDgp(grid[j], 1, mean_hat, sd_hat) / (0.5 * n^2)) } sd_dgp = sd_dgp :* sqrt(abs((invsym(S) * S2 * invsym(S))[v+1, v+1])) } // bandwidth S = optimize_init() optimize_init_evaluator(S, &lpdensity_optFunc()) optimize_init_params(S, 1e-8) optimize_init_argument(S, 1, bias_dgp) optimize_init_argument(S, 2, sd_dgp) optimize_init_argument(S, 3, p) optimize_init_argument(S, 4, v) optimize_init_argument(S, 5, "imse-rot") optimize_init_which(S, "min") optimize_init_evaluatortype(S, "gf0") optimize_init_tracelevel(S, "none") optimize_init_conv_maxiter(S, 20) (void) optimize(S) h = optimize_result_params(S) if (regularize == 1) { for (j=1; j<=ng; j++) { if (nLocalMin > 0) { h = max((h, sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))])) } if (nUniqueMin > 0) { h = max((h, sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))])) } } } h = min(( h, max(( abs(max(dataUnique)-min(grid)), abs(min(dataUnique)-max(grid)) )) )) if (stdVar) { h = h:* scale_temp } return(h) } mata mosave lpdensity_bwIROT(), replace end ******************************************************************************** * mse-dpi bandwidth ******************************************************************************** capture mata mata drop lpdensity_bwMSE() mata real matrix lpdensity_bwMSE(real colvector data, real colvector grid, real scalar p, real scalar v, string scalar kernel, real colvector Cweights, real colvector Pweights, real scalar massPoints, real scalar stdVar, real scalar regularize, real scalar nLocalMin, real scalar nUniqueMin){ ii = order(data, 1) data = data[ii, 1] Cweights = Cweights[ii, 1] Pweights = Pweights[ii, 1] n = length(data) ng = length(grid) dataUnique = lpdensity_unique(data) freqUnique = dataUnique[., 2] indexUnique = dataUnique[., 3] dataUnique = dataUnique[., 1] nUnique = length(dataUnique) if (stdVar) { center_temp = mean(data) scale_temp = sqrt(variance(data)) data = (data :- center_temp) :/ scale_temp dataUnique = (dataUnique :- center_temp) :/ scale_temp grid = (grid :- center_temp) :/ scale_temp } if (massPoints) { Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique) } else{ Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights) } weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n Cweights = Cweights :/ sum(Cweights) :* n Pweights = Pweights :/ sum(Pweights) :* n weights_normalUnique = runningsum(weights_normal)[indexUnique] if (nUnique > 1) { weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) } CweightsUnique = runningsum(Cweights)[indexUnique] if (nUnique > 1) { CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) } PweightsUnique = runningsum(Pweights)[indexUnique] if (nUnique > 1) { PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) } // obtain preliminary bandwidth for estimating densities // this is used for constructing preasymptotic matrices h1 = lpdensity_bwIROT(data, grid, 2, 1, kernel, Cweights, Pweights, 1, 1, 1, 20+2+1, 20+2+1) // obtain preliminary bandwidth for estimating F_p+1 // this is used for constructing F_p+1 hp1 = lpdensity_bwIROT(data, grid, p+2, p+1, kernel, Cweights, Pweights, 1, 1, 1, 20+p+2+1, 20+p+2+1) // obtain preliminary bandwidth for estimating F_p+2 // this is used for constructing F_p+2 hp2 = lpdensity_bwIROT(data, grid, p+3, p+2, kernel, Cweights, Pweights, 1, 1, 1, 20+p+3+1, 20+p+3+1) dgp_hat = J(ng, 2, .) // Fp+1 and Fp+2 with normalization constants const_hat = J(ng, 3, .) h = J(ng, 1, .) for (j=1; j<=ng; j++) { // estimate F_p+2 index_temp = (abs(data :- grid[j]) :<= hp2) Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp2 Xh_p_temp = J(sum(index_temp), p+4, 1) Xh_p_temp[., 2] = Xh_temp for (jj=3; jj<=p+4; jj++) { Xh_p_temp[., jj] = Xh_temp :^ (jj-1) } if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ hp2 if (kernel == "uniform") Kh_temp = 0.5 / hp2 if (kernel == "epanechnikov") Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp2 Kh_temp = select(Pweights, index_temp) :* Kh_temp Y_temp = select(Fn, index_temp) dgp_hat[j, 2] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+3, 1] / (hp2^(p+2)) // estimate F_p+1 index_temp = (abs(data :- grid[j]) :<= hp1) Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp1 Xh_p_temp = J(sum(index_temp), p+3, 1) Xh_p_temp[., 2] = Xh_temp for (jj=3; jj<=p+3; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) } if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ hp1 if (kernel == "uniform") Kh_temp = 0.5 / hp1 if (kernel == "epanechnikov") Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp1 Kh_temp = select(Pweights, index_temp) :* Kh_temp Y_temp = select(Fn, index_temp) dgp_hat[j, 1] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+2, 1] / (hp1^(p+1)) // prepare for estimating matrices index_temp = (abs(data :- grid[j]) :<= h1) Xh_temp = (select(data, index_temp) :- grid[j]) :/ h1 if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ h1 if (kernel == "uniform") Kh_temp = 0.5 / h1 if (kernel == "epanechnikov") Kh_temp = 0.75 :* (1 :- Xh_temp:^2) :/ h1 Pweights_temp = select(Pweights, index_temp) // estimate Cp matrix Xh_p_temp = J(sum(index_temp), (2*p+1) - (p+1) + 1, 1) for (jj=1; jj<=(2*p+1)-(p+1)+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(p+jj) } C_p_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n // estimate Cp+1 matrix Xh_p_temp = J(sum(index_temp), (2*p+2) - (p+2) + 1, 1) for (jj=1; jj<=(2*p+2)-(p+2)+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(p+jj+1) } C_p1_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n // estimate S matirx Xh_p_temp = J(sum(index_temp), p - 0 + 1, 1) for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) } S_hat = cross(Xh_p_temp, Kh_temp :* Pweights_temp, Xh_p_temp) :/ n S_hat_inv = invsym(S_hat) // estimate G matrix if (v == 0) { G_hat = cross(Xh_p_temp, Kh_temp:^2 :* Pweights_temp, Xh_p_temp) :/ n } else { if (massPoints) { Y_temp = Fn[indexUnique] Xh_temp = (dataUnique :- grid[j]) :/ h1 Xh_p_temp = J(nUnique, p - 0 + 1, 1) Xh_p_Kh_Pweights_temp = Xh_p_temp if (kernel == "triangular") Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp[indexUnique] if (kernel == "uniform") Kh_temp = (0.5 / h1) :* index_temp[indexUnique] if (kernel == "epanechnikov") Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp[indexUnique] for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique } F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1;jj++) { G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal } G_hat = cross(G, G) :/ n } else { Y_temp = Fn Xh_temp = (data :- grid[j]) :/ h1 Xh_p_temp = J(n, p - 0 + 1, 1) Xh_p_Kh_Pweights_temp = Xh_p_temp if (kernel == "triangular") Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp if (kernel == "uniform") Kh_temp = (0.5 / h1) :* index_temp if (kernel == "epanechnikov") Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights } F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1;jj++) { G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal } G_hat = cross(G, G) :/ n } } // now get all constants const_hat[j, 1] = factorial(v) * (S_hat_inv * C_p_hat )[v+1] const_hat[j, 2] = factorial(v) * (S_hat_inv * C_p1_hat)[v+1] if (v > 0) { const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (n * h1)) } else { temp_ii = min(( max(( mean(data :<= grid[j]), 1/n )), 1 - 1/n )) const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (0.5 * n^2) * h1 * temp_ii * (1 - temp_ii)) } // now optimal bandwidth S = optimize_init() optimize_init_evaluator(S, &lpdensity_optFunc()) optimize_init_params(S, 1e-8) optimize_init_argument(S, 1, dgp_hat[j, .]) optimize_init_argument(S, 2, const_hat[j, .]) optimize_init_argument(S, 3, p) optimize_init_argument(S, 4, v) optimize_init_argument(S, 5, "mse-dpi") optimize_init_which(S, "min") optimize_init_evaluatortype(S, "gf0") optimize_init_tracelevel(S, "none") optimize_init_conv_maxiter(S, 20) (void) optimize(S) h[j] = optimize_result_params(S) if (regularize == 1) { if (nLocalMin > 0) { h[j] = max((h[j], sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))])) } if (nUniqueMin > 0) { h[j] = max((h[j], sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))])) } } h[j] = min(( h[j], max( abs(dataUnique :- grid[j]) ) )) } if (stdVar) { h = h :* scale_temp } return(h) } mata mosave lpdensity_bwMSE(), replace end ******************************************************************************** * imse-dpi bandwidth ******************************************************************************** capture mata mata drop lpdensity_bwIMSE() mata real scalar lpdensity_bwIMSE(real colvector data, real colvector grid, real scalar p, real scalar v, string scalar kernel, real colvector Cweights, real colvector Pweights, real scalar massPoints, real scalar stdVar, real scalar regularize, real scalar nLocalMin, real scalar nUniqueMin){ ii = order(data, 1) data = data[ii, 1] Cweights = Cweights[ii, 1] Pweights = Pweights[ii, 1] n = length(data) ng = length(grid) dataUnique = lpdensity_unique(data) freqUnique = dataUnique[., 2] indexUnique = dataUnique[., 3] dataUnique = dataUnique[., 1] nUnique = length(dataUnique) if (stdVar) { center_temp = mean(data) scale_temp = sqrt(variance(data)) data = (data :- center_temp) :/ scale_temp dataUnique = (dataUnique :- center_temp) :/ scale_temp grid = (grid :- center_temp) :/ scale_temp } if (massPoints) { Fn = lpdensity_rep((runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights))[indexUnique], freqUnique) } else{ Fn = runningsum(Cweights :* Pweights) :/ sum(Cweights :* Pweights) } weights_normal = Cweights :* Pweights :/ sum(Cweights :* Pweights) :* n Cweights = Cweights :/ sum(Cweights) :* n Pweights = Pweights :/ sum(Pweights) :* n weights_normalUnique = runningsum(weights_normal)[indexUnique] if (nUnique > 1) { weights_normalUnique = weights_normalUnique :- (0\ weights_normalUnique[1..(nUnique-1)]) } CweightsUnique = runningsum(Cweights)[indexUnique] if (nUnique > 1) { CweightsUnique = CweightsUnique :- (0\ CweightsUnique[1..(nUnique-1)]) } PweightsUnique = runningsum(Pweights)[indexUnique] if (nUnique > 1) { PweightsUnique = PweightsUnique :- (0\ PweightsUnique[1..(nUnique-1)]) } // obtain preliminary bandwidth for estimating densities // this is used for constructing preasymptotic matrices h1 = lpdensity_bwIROT(data, grid, 2, 1, kernel, Cweights, Pweights, 1, 1, 1, 20+2+1, 20+2+1) // obtain preliminary bandwidth for estimating F_p+1 // this is used for constructing F_p+1 hp1 = lpdensity_bwIROT(data, grid, p+2, p+1, kernel, Cweights, Pweights, 1, 1, 1, 20+p+2+1, 20+p+2+1) // obtain preliminary bandwidth for estimating F_p+2 // this is used for constructing F_p+2 hp2 = lpdensity_bwIROT(data, grid, p+3, p+2, kernel, Cweights, Pweights, 1, 1, 1, 20+p+3+1, 20+p+3+1) dgp_hat = J(ng, 2, .) // Fp+1 and Fp+2 with normalization constants const_hat = J(ng, 3, .) for (j=1; j<=ng; j++) { // estimate F_p+2 index_temp = (abs(data :- grid[j]) :<= hp2) Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp2 Xh_p_temp = J(sum(index_temp), p+4, 1) Xh_p_temp[., 2] = Xh_temp for (jj=3; jj<=p+4; jj++) { Xh_p_temp[., jj] = Xh_temp :^ (jj-1) } if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ hp2 if (kernel == "uniform") Kh_temp = 0.5 / hp2 if (kernel == "epanechnikov") Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp2 Kh_temp = select(Pweights, index_temp) :* Kh_temp Y_temp = select(Fn, index_temp) dgp_hat[j, 2] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+3, 1] / (hp2^(p+2)) // estimate F_p+1 index_temp = (abs(data :- grid[j]) :<= hp1) Xh_temp = (select(data, index_temp) :- grid[j]) :/ hp1 Xh_p_temp = J(sum(index_temp), p+3, 1) Xh_p_temp[., 2] = Xh_temp for (jj=3; jj<=p+3; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) } if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ hp1 if (kernel == "uniform") Kh_temp = 0.5 / hp1 if (kernel == "epanechnikov") Kh_temp = (1 :- (Xh_temp):^2) :* 0.75 :/ hp1 Kh_temp = select(Pweights, index_temp) :* Kh_temp Y_temp = select(Fn, index_temp) dgp_hat[j, 1] = (invsym(cross(Xh_p_temp, Kh_temp, Xh_p_temp)) * cross(Xh_p_temp, Kh_temp, Y_temp))[p+2, 1] / (hp1^(p+1)) // prepare for estimating matrices index_temp = (abs(data :- grid[j]) :<= h1) Xh_temp = (select(data, index_temp) :- grid[j]) :/ h1 if (kernel == "triangular") Kh_temp = (1 :- abs(Xh_temp)) :/ h1 if (kernel == "uniform") Kh_temp = 0.5 / h1 if (kernel == "epanechnikov") Kh_temp = 0.75 :* (1 :- Xh_temp:^2) :/ h1 Pweights_temp = select(Pweights, index_temp) // estimate Cp matrix Xh_p_temp = J(sum(index_temp), (2*p+1) - (p+1) + 1, 1) for (jj=1; jj<=(2*p+1)-(p+1)+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(p+jj) } C_p_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n // estimate Cp+1 matrix Xh_p_temp = J(sum(index_temp), (2*p+2) - (p+2) + 1, 1) for (jj=1; jj<=(2*p+2)-(p+2)+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(p+jj+1) } C_p1_hat = cross(Xh_p_temp, Kh_temp, Pweights_temp) :/ n // estimate S matirx Xh_p_temp = J(sum(index_temp), p - 0 + 1, 1) for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) } S_hat = cross(Xh_p_temp, Kh_temp :* Pweights_temp, Xh_p_temp) :/ n S_hat_inv = invsym(S_hat) // estimate G matrix if (v == 0) { G_hat = cross(Xh_p_temp, Kh_temp:^2 :* Pweights_temp, Xh_p_temp) :/ n } else { if (massPoints) { Y_temp = Fn[indexUnique] Xh_temp = (dataUnique :- grid[j]) :/ h1 Xh_p_temp = J(nUnique, p - 0 + 1, 1) Xh_p_Kh_Pweights_temp = Xh_p_temp if (kernel == "triangular") Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp[indexUnique] if (kernel == "uniform") Kh_temp = (0.5 / h1) :* index_temp[indexUnique] if (kernel == "epanechnikov") Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp[indexUnique] for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* PweightsUnique } F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1;jj++) { G[., jj] = (lpdensity_rep((runningsum(Xh_p_Kh_Pweights_temp[nUnique..1, jj]) :/ n)[nUnique..1], freqUnique) :- F_Xh_p_Kh_temp[jj]) :* weights_normal } G_hat = cross(G, G) :/ n } else { Y_temp = Fn Xh_temp = (data :- grid[j]) :/ h1 Xh_p_temp = J(n, p - 0 + 1, 1) Xh_p_Kh_Pweights_temp = Xh_p_temp if (kernel == "triangular") Kh_temp = ((1 :- abs(Xh_temp)) :/ h1) :* index_temp if (kernel == "uniform") Kh_temp = (0.5 / h1) :* index_temp if (kernel == "epanechnikov") Kh_temp = (0.75 :* (1 :- Xh_temp:^2) :/ h1) :* index_temp for (jj=1; jj<=p-0+1; jj++) { Xh_p_temp[., jj] = Xh_temp:^(jj-1) Xh_p_Kh_Pweights_temp[., jj] = Xh_p_temp[., jj] :* Kh_temp :* Pweights } F_Xh_p_Kh_temp = cross(Xh_p_Kh_Pweights_temp, Y_temp) :/ n G = J(n, p - 0 + 1, .) for (jj=1; jj<=p-0+1;jj++) { G[., jj] = ((runningsum(Xh_p_Kh_Pweights_temp[n..1, jj]) :/ n)[n..1] :- F_Xh_p_Kh_temp[jj]) :* weights_normal } G_hat = cross(G, G) :/ n } } // now get all constants const_hat[j, 1] = factorial(v) * (S_hat_inv * C_p_hat )[v+1] const_hat[j, 2] = factorial(v) * (S_hat_inv * C_p1_hat)[v+1] if (v > 0) { const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (n * h1)) } else { temp_ii = min(( max(( mean(data :<= grid[j]), 1/n )), 1 - 1/n )) const_hat[j, 3] = factorial(v) * sqrt(abs((S_hat_inv * G_hat * S_hat_inv)[v+1,v+1]) / (0.5 * n^2) * h1 * temp_ii * (1 - temp_ii)) } } // bandwidth S = optimize_init() optimize_init_evaluator(S, &lpdensity_optFunc()) optimize_init_params(S, 1e-8) optimize_init_argument(S, 1, dgp_hat) optimize_init_argument(S, 2, const_hat) optimize_init_argument(S, 3, p) optimize_init_argument(S, 4, v) optimize_init_argument(S, 5, "imse-dpi") optimize_init_which(S, "min") optimize_init_evaluatortype(S, "gf0") optimize_init_tracelevel(S, "none") optimize_init_conv_maxiter(S, 20) (void) optimize(S) h = optimize_result_params(S) if (regularize == 1) { for (j=1; j<=ng; j++) { if (nLocalMin > 0) { h = max((h, sort(abs(data:-grid[j]), 1)[min((n, nLocalMin))])) } if (nUniqueMin > 0) { h = max((h, sort(abs(dataUnique:-grid[j]), 1)[min((nUnique, nUniqueMin))])) } } } h = min(( h, max(( abs(max(dataUnique)-min(grid)), abs(min(dataUnique)-max(grid)) )) )) if (stdVar) { h = h :* scale_temp } return(h) } mata mosave lpdensity_bwIMSE(), replace end ******************************************************************************** * optimization criterion function ******************************************************************************** capture mata mata drop lpdensity_optFunc() mata void lpdensity_optFunc(todo, h, mat1, mat2, p, v, bwselect, fn, gn, hn) { if (bwselect == "imse-dpi") { if (v > 0) { fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :* mat2[., 1] :+ h :* mat1[., 2] :* mat2[., 2]):^2 + mat2[., 3]:^2 :/ h^(2*v - 1) ) } else { fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :* mat2[., 1] :+ h :* mat1[., 2] :* mat2[., 2]):^2 + mat2[., 3]:^2 :/ h ) } } else if (bwselect == "imse-rot") { if (v > 0) { fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :+ h :* mat1[., 2]):^2 + mat2[., 1]:^2 :/ h^(2*v - 1) ) } else { fn = sum( h^(2*p+2-2*v) :* (mat1[., 1] :+ h :* mat1[., 2]):^2 + mat2[., 1]:^2 :/ h ) } } else if (bwselect == "mse-dpi") { if (v > 0) { fn = h^(2*p+2-2*v) * (mat1[1] * mat2[1] + h * mat1[2] * mat2[2])^2 + mat2[3]^2 / h^(2*v - 1) } else { fn = h^(2*p+2-2*v) * (mat1[1] * mat2[1] + h * mat1[2] * mat2[2])^2 + mat2[3]^2 / h } } else { if (v > 0) { fn = h^(2*p+2-2*v) * (mat1[1] + h * mat1[2])^2 + mat2[1]^2 / h^(2*v - 1) } else { fn = h^(2*p+2-2*v) * (mat1[1] + h * mat1[2])^2 + mat2[1]^2 / h } } } mata mosave lpdensity_optFunc(), replace end