//////////////////////////////////////////////////////////////////////////////// // STATA FOR Beresteanu, A. & Sasaki, Y. (2021): Quantile Regression with // Interval Data. Econometric Reviews (Special Issue Honoring Cheng // Hsiao), 40 (6): 562-583. // // Use it when you have interval-valued data (consisting of lower and upper // bounds) and want to compute set percentiles with their confidence sets. //////////////////////////////////////////////////////////////////////////////// program define itvalpctile, rclass version 14.2 syntax varlist(min=2 max=2 numeric) [if] [in] [,cover(real 0.95) pl(real 10) ph(real 90) np(real 9) conditional(varname numeric) location(real 999999)] marksample touse gettoken depvar indepvars : varlist _fv_check_depvar `depvar' fvexpand `indepvars' local cnames `r(varlist)' tempname N percent lowerEST upperEST lowerCS upperCS if "`conditional'" == "" { mata: estimate_setp("`depvar'", "`cnames'", `pl', `ph', `np', "`touse'", `cover', "`N'", /// "`percent'", "`lowerEST'", "`upperEST'", "`lowerCS'", "`upperCS'") } if "`conditional'" != "" { tempvar x gen `x' = `conditional' mata: estimate_conditional_setp("`depvar'", "`cnames'", "`x'", `pl', `ph', `np', "`touse'", `cover', `location', "`N'", /// "`percent'", "`lowerEST'", "`upperEST'", "`lowerCS'", "`upperCS'") } return scalar N = `N' return local cmd "itvalpctile" return matrix CSupper = `upperCS' return matrix CSlower = `lowerCS' return matrix upper = `upperEST' return matrix lower = `lowerEST' return matrix percent = `percent' end mata: //////////////////////////////////////////////////////////////////////////////// // Kernel Function void kernel(u, kout){ kout = 1 :* ( -1 :< u ) :* ( u :< 1 ) } //////////////////////////////////////////////////////////////////////////////// // Unconditional void estimate_setp( string scalar yv, string scalar xv, real scalar q_low, real scalar q_high, real scalar q_num, string scalar touse, real scalar cover, string scalar nname, string scalar pname, string scalar lename, string scalar uename, string scalar lcname, string scalar ucname) { printf("\n{hline 78}\n") printf("Executing: Beresteanu, A. & Sasaki, Y. (2020): Quantile Regression with\n") printf(" Interval Data. Econometric Reviews (Special Issue Honoring Cheng\n") printf(" Hsiao), 40 (6): 562-583.\n") printf("{hline 78}\n") yl = st_data(., yv, touse) yh = st_data(., xv, touse) n = rows(yl) original_n = n // List of quantiles at which to evaluate causal effects qlist = (q_high:/100 - q_low:/100) :* (0..(q_num-1)) :/ (q_num-1) :+ q_low:/100 //////////////////////////////////////////////////////////////////////////// // Estimate interval quantiles lower = sort(yl,1)[trunc(n:*qlist)] upper = sort(yh,1)[trunc(n:*qlist)] //////////////////////////////////////////////////////////////////////////// // Confidence Sets real matrix kout1, kout2 conf_lower = lower conf_upper = upper h1 = 1.06 * variance(yl)^0.5 / n^0.2 h2 = 1.06 * variance(yh)^0.5 / n^0.2 for( idx = 1 ; idx <= length(qlist) ; idx++ ){ kernel( (yl :- lower[idx]) :/ h1, kout1 ) kernel( (yh :- upper[idx]) :/ h2, kout2 ) fa = sum( kout1 ) / (n * h1) fb = sum( kout2 ) / (n * h2) Fab = mean( yl :<= lower[idx] :& yh :<= upper[idx] ) Sigma11 = qlist[idx] * (1-qlist[idx]) / fa^2 Sigma12 = (Fab - qlist[idx]^2) / (fa * fb) Sigma22 = qlist[idx] * (1-qlist[idx]) / fb^2 scale = invnormal(1 - (1 - cover) / 2) conf_lower[idx] = lower[idx] - scale * ( Sigma11 / n )^0.5 conf_upper[idx] = upper[idx] + scale * ( Sigma22 / n )^0.5 } //////////////////////////////////////////////////////////////////////////// // Console output printf("\n{hline 78}\n") printf( "Number of observations: n = %f\n", original_n) printf( "{hline 78}\n") printf( "Percent [ Percentile Itervals ] [ %2.0f%% Confidence Sets ]\n",100*cover) printf( "{hline 71}\n") for( idx = 1 ; idx <= length(qlist) ; idx++ ){ printf(" %3.0f [%10.4f %10.4f] [%10.4f %10.4f]\n",qlist[idx]*100,lower[idx],upper[idx],conf_lower[idx],conf_upper[idx]) } printf( "{hline 78}\n") st_numscalar(nname, n) st_matrix(pname, qlist':*100) st_matrix(lename, lower) st_matrix(uename, upper) st_matrix(lcname, conf_lower) st_matrix(ucname, conf_upper) } //////////////////////////////////////////////////////////////////////////////// // Conditional void estimate_conditional_setp( string scalar yv, string scalar xv, string scalar xx, real scalar q_low, real scalar q_high, real scalar q_num, string scalar touse, real scalar cover, real scalar location, string scalar nname, string scalar pname, string scalar lename, string scalar uename, string scalar lcname, string scalar ucname) { printf("\n{hline 78}\n") printf("Executing: Beresteanu, A. & Sasaki, Y. (2020): Quantile Regression with\n") printf(" Interval Data. Econometric Reviews (Special Issue Honoring Cheng\n") printf(" Hsiao), 40 (6): 562-583.\n") printf("{hline 78}\n") yl = st_data(., yv, touse) yh = st_data(., xv, touse) x = st_data(., xx, touse) n = rows(yl) original_n = n // List of quantiles at which to evaluate causal effects qlist = (q_high:/100 - q_low:/100) :* (0..(q_num-1)) :/ (q_num-1) :+ q_low:/100 //////////////////////////////////////////////////////////////////////////// // Find local observations if( location == 999999 ){ location = mean(x) } h = 1.06 * variance(x)^0.5 / n^0.2 yl = select(yl, abs( x :- location ) :<= h) yh = select(yh, abs( x :- location ) :<= h) n = rows(yl) //////////////////////////////////////////////////////////////////////////// // Estimate interval quantiles lower = sort(yl,1)[trunc(n:*qlist)] upper = sort(yh,1)[trunc(n:*qlist)] //////////////////////////////////////////////////////////////////////////// // Confidence Sets real matrix kout1, kout2 conf_lower = lower conf_upper = upper h1 = 1.06 * variance(yl)^0.5 / n^0.2 h2 = 1.06 * variance(yh)^0.5 / n^0.2 for( idx = 1 ; idx <= length(qlist) ; idx++ ){ kernel( (yl :- lower[idx]) :/ h1, kout1 ) kernel( (yh :- upper[idx]) :/ h2, kout2 ) fa = sum( kout1 ) / (n * h1) fb = sum( kout2 ) / (n * h2) Fab = mean( yl :<= lower[idx] :& yh :<= upper[idx] ) Sigma11 = qlist[idx] * (1-qlist[idx]) / fa^2 Sigma12 = (Fab - qlist[idx]^2) / (fa * fb) Sigma22 = qlist[idx] * (1-qlist[idx]) / fb^2 scale = invnormal(1 - (1 - cover) / 2) conf_lower[idx] = lower[idx] - scale * ( Sigma11 / n )^0.5 conf_upper[idx] = upper[idx] + scale * ( Sigma22 / n )^0.5 } //////////////////////////////////////////////////////////////////////////// // Console output printf("\n{hline 78}\n") printf( "Number of observations: n = %f\n", original_n) printf( "Localition of the conditioning variable: x = %f\n", location) printf( "{hline 78}\n") printf( "Percent [ Percentile Itervals ] [ %2.0f%% Confidence Sets ]\n",100*cover) printf( "{hline 71}\n") for( idx = 1 ; idx <= length(qlist) ; idx++ ){ printf(" %3.0f [%10.4f %10.4f] [%10.4f %10.4f]\n",qlist[idx]*100,lower[idx],upper[idx],conf_lower[idx],conf_upper[idx]) } printf( "{hline 78}\n") st_numscalar(nname, n) st_matrix(pname, qlist':*100) st_matrix(lename, lower) st_matrix(uename, upper) st_matrix(lcname, conf_lower) st_matrix(ucname, conf_upper) } end ////////////////////////////////////////////////////////////////////////////////