//////////////////////////////////////////////////////////////////////////////// // STATA FOR Chiang, H.D. & Sasaki, Y. (2019): Causal Inference by Quantile // Regression Kink Designs. Journal of Econometrics 210 (2), 405-433. // // Use it when you consider a regression kink design and you are interested in // analyzing heterogeneous causal effects (e.g., heterogeneous effects of // unemployment insurance on unemployment duration). //////////////////////////////////////////////////////////////////////////////// program define qrkd, rclass version 14.2 syntax varlist(numeric) [if] [in] [, k(real 0) bpl(real 0) bpr(real 1) cover(real 0.95) ql(real 0.25) qh(real 0.75) qn(real 3) bw(real -1)] marksample touse gettoken depvar indepvars : varlist _fv_check_depvar depvar' fvexpand indepvars' local cnames r(varlist)' tempname b V N cb h q bl bu mata: estimate_qrkd("depvar'", "cnames'", /// k', bpr', bpl', /// ql', qh', qn', /// bw', "touse'", /// "b'", "V'", "N'", /// cover', "cb'", "h'", /// "q'", "bl'", "bu'") matrix colnames b' = QRKD matrix colnames V' = QRKD matrix rownames V' = QRKD return scalar cover= cover' return scalar k = k' return scalar h = h' return scalar N = N' return matrix V = V' return matrix CBupper = bu' return matrix CBlower = bl' return matrix b = b' return matrix q = q' return local cmd "qrkd" end mata: //////////////////////////////////////////////////////////////////////////////// // Kernel Function void kernel(u, kout){ kout = (70/81) :* (1 :- (u:^2):^(3/2) ):^3 :* ( -1 :< u ) :* ( u :< 1 ) } //////////////////////////////////////////////////////////////////////////////// // Smoothed Check Function void check(u, tau, smooth, checkout){ ch1 = (u :- smooth) :* (tau :- (u :- smooth :< 0) ) ch2 = (u :- smooth:/2) :* (tau :- (u :- smooth:/2 :< 0) ) ch3 = (u) :* (tau :- (u :< 0) ) ch4 = (u :+ smooth:/2) :* (tau :- (u :+ smooth:/2 :< 0) ) ch5 = (u :+ smooth) :* (tau :- (u :+ smooth :< 0) ) checkout = ch1 :/ 5 + ch2 :/ 5 + ch3 :/ 5 + ch4 :/ 5 + ch5 :/ 5 } //////////////////////////////////////////////////////////////////////////////// // Function for the Criterion void qcrit(todo, para, y, x, h, tau, crit, g, H){ alph = para[1] beta1plus = para[2] beta1minus = para[3] beta2plus = para[4] beta2minus = para[5] real matrix checkout real matrix kout smooth = variance(x)^0.5 / 1000 check(y :- alph :- beta1plus :* x :* (x :> 0) :- beta1minus :* x :* (x :< 0) :- beta2plus :* x:^2 :* (x :> 0) :/ 2 :- beta2minus :* x:^2 :* (x :< 0) :/ 2, tau, smooth, checkout) kernel(x:/h, kout) crit = sum( kout :* checkout ) + rows(y) * 0.001 * sum( beta2plus^2 + beta2minus^2 ) } //////////////////////////////////////////////////////////////////////////////// // Estimation of QRKD void estimate_qrkd( string scalar yv, string scalar xv, real scalar cut, real scalar bpright, real scalar bpleft, real scalar q_low, real scalar q_high, real scalar q_num, real scalar b_w, string scalar touse, string scalar bname, string scalar Vname, string scalar nname, real scalar cover, string scalar cbname, string scalar hname, string scalar qname, string scalar blname, string scalar buname) { printf("\n{hline 78}\n") printf("Executing: Chiang, H.D. & Sasaki, Y. (2019): Causal Inference by Quantile\n") printf(" Regression Kink Designs. Journal of Econometrics 210 (2), 405-433.\n") printf("{hline 78}\n") real vector y, x, beta1plus_hat, beta1minus_hat, qrkd_hat real scalar cutoff, n y = st_data(., yv, touse) x = st_data(., xv, touse) :- cut n = rows(y) // List of quantiles at which to evaluate causal effects qlist = (q_high - q_low) :* (0..(q_num-1)) :/ (q_num-1) :+ q_low //////////////////////////////////////////////////////////////////////////// // Estimate QRKD init_alpha = J(length(qlist),1,0) alph_hat = J(length(qlist),1,0) beta1plus_hat = J(length(qlist),1,0) beta1minus_hat = J(length(qlist),1,0) qrkd_hat = J(length(qlist),1,0) for( idx = 1 ; idx <= length(qlist) ; idx++ ){ printf("Iteration %f/%f: Estimating QRKD at the %f-th quantile\n", idx, length(qlist)+1, qlist[idx]) h = b_w if( b_w <= 0 ){ h = 10 * variance(x)^0.5 * variance(y)^0.5 / n^0.2 } tau = qlist[idx] //////////////////////////////////////////////////////////////////////// // Set initial parameter values xleft = select(x, x :<= 0) xleft = J(rows(xleft),1,1), xleft xright = select(x, x :> 0) xright = J(rows(xright),1,1), xright yleft = select(y, x :<= 0) yright = select(y, x :> 0) yclose = select(y, -variance(x):^0.5:/2 :<= x :& x :<= variance(x):^0.5:/2) init_alpha[idx] = sort(yclose,1)[trunc(rows(yclose)*tau)] beta1plus = (luinv(xright'*xright)*xright'*yright)[2,1] beta1minus = (luinv(xleft'*xleft)*xleft'*yleft)[2,1] beta2plus = 0 beta2minus = 0 initpara = (init_alpha[idx], beta1plus, beta1minus, beta2plus, beta2minus) //////////////////////////////////////////////////////////////////////// // Optimization routine S = optimize_init() optimize_init_evaluator(S,&qcrit()) optimize_init_which(S,"min") optimize_init_evaluatortype(S, "d0") optimize_init_technique(S,"nr") optimize_init_singularHmethod(S,"hybrid") optimize_init_argument(S,1,y) optimize_init_argument(S,2,x) optimize_init_argument(S,3,h) optimize_init_argument(S,4,tau) optimize_init_tracelevel(S, "none") optimize_init_conv_warning(S, "off") optimize_init_params(S, initpara) est=optimize(S) alph_hat[idx] = est[1] beta1plus_hat[idx] = est[2] beta1minus_hat[idx] = est[3] qrkd_hat[idx] = (beta1plus_hat[idx] - beta1minus_hat[idx]) / (bpright - bpleft) } //////////////////////////////////////////////////////////////////////////// // Compute K_{in\tau} and z_{in\tau} real matrix Kintau, zintau kernel(x:/h, Kintau) Kintau = Kintau, Kintau, Kintau, Kintau, Kintau zintau = J(n,1,1), x :/ h :* (x :> 0), x :/ h :* (x :<= 0), (x :/ h):^2 :* (x :> 0), (x :/ h):^2 :* (x :<= 0) //////////////////////////////////////////////////////////////////////////// // Compute the N matrix Nmatrix = J(5,5,0) ulist = (-100..100) :/ 100 ulistinterval = ulist[1] - ulist[2] real matrix kernelulist kernel(ulist, kernelulist) for( idx = 1 ; idx <= length(ulist) ; idx++ ){ u = ulist[idx] ubar = 1 \ u * (u > 0) \ u * (u <= 0) \ u^2 * (u > 0) \ u^2 * (u <= 0) Nmatrix = Nmatrix + ( (ubar * ubar') :* (ulistinterval * kernelulist[idx]) ) } //////////////////////////////////////////////////////////////////////////// // Compute other auxiliary objects for variance estimation i2minusi3 = (0 \ 1 \ -1 \ 0 \ 0) corr_xy = mean( x :* y ) - mean(x) * mean(y) corr_xy = corr_xy / (variance(x) * variance(y))^0.5 hx = variance(x)^0.5 / n^(1/5) hxx = 5 * variance(x)^(1.5/2) * variance(y)^(1/2) / n^(1/6) // 10 * variance(x)^0.5 / n^(1/6) real matrix kxlist, kxxlist kernel(x:/hx, kxlist) kernel(x:/hxx, kxxlist) fx = mean(kxlist)/hx hy = hxx / variance(x)^(1/2) //10 * variance(y)^0.5 / n^(1/6) fyx = J(length(qlist),1,1) for( idx = 1 ; idx <= length(qlist) ; idx++ ){ real matrix kylist kernel((y :- alph_hat[idx]):/hy, kylist) fyx[idx] = ( mean(kxxlist :* kylist) / (hxx * hy) ) / ( mean(kxxlist) / hxx ) } //////////////////////////////////////////////////////////////////////////// // Pivotal approach to variance estimation printf("Iteration %f/%f: Variance Estimation\n", length(qlist)+1, length(qlist)+1) num_bootstrap = 4000 Ylist = J(length(qlist),num_bootstrap,0) for( idx = 1 ; idx <= num_bootstrap ; idx++ ){ unif = uniform(n,1) unif = unif, unif, unif, unif, unif for( jdx = 1 ; jdx <= length(qlist) ; jdx++ ){ numerator = i2minusi3' * luinv(Nmatrix) * ( zintau :* Kintau :* (qlist[jdx] :- 1:*(unif :<= qlist[jdx])) )' * J(n,1,1) denominator = (bpright - bpleft) * (n*h)^0.5 * fx * fyx[jdx] Ylist[jdx,idx] = numerator / denominator } } sigma = ( diagonal(Ylist * Ylist' :/ num_bootstrap - (Ylist :/ num_bootstrap) * (Ylist :/ num_bootstrap)') :/ (n*h^3) ):^0.5 //////////////////////////////////////////////////////////////////////////// // Standardized pivotal approach Ystdlist = J(length(qlist),num_bootstrap,0) for( idx = 1 ; idx <= num_bootstrap ; idx++ ){ unif = uniform(n,1) unif = unif, unif, unif, unif, unif for( jdx = 1 ; jdx <= length(qlist) ; jdx++ ){ numerator = i2minusi3' * luinv(Nmatrix) * ( zintau :* Kintau :* (qlist[jdx] :- 1:*(unif :<= qlist[jdx])) )' * J(n,1,1) denominator = (bpright - bpleft) * (n*h)^0.5 * fx * fyx[jdx] Ystdlist[jdx,idx] = numerator / denominator / sigma[jdx] } } //////////////////////////////////////////////////////////////////////////// // Hypothesis testing WS = (n*h^3)^0.5 * max(abs(qrkd_hat)) WH = (n*h^3)^0.5 * max(abs( qrkd_hat :- mean(qrkd_hat) )) WSboot = J(num_bootstrap,1,0) WHboot = J(num_bootstrap,1,0) for( idx = 1 ; idx <= num_bootstrap ; idx++ ){ WSboot[idx] = max(abs(Ylist[.,idx])) WHboot[idx] = max(abs( Ylist[.,idx] :- mean(Ylist[.,idx]) )) } pvalueWS = mean( 1 * (WSboot :> WS) ) pvalueWH = mean( 1 * (WHboot :> WH) ) //////////////////////////////////////////////////////////////////////////// // Uniform confidence band WSbootStd = J(num_bootstrap,1,0) for( idx = 1 ; idx <= num_bootstrap ; idx++ ){ WSbootStd[idx] = max(abs(Ystdlist[.,idx])) } halfbandlength = ( sort(WSbootStd,1)[trunc(cover * num_bootstrap)] / (n*h^3)^0.5 ) :* sigma //////////////////////////////////////////////////////////////////////////// // Set b, V, n and cb b = qrkd_hat' V = Ylist * Ylist' :/ num_bootstrap - (Ylist :/ num_bootstrap) * (Ylist :/ num_bootstrap)' V = V :/ (n*h^3) bu = qrkd_hat' :+ halfbandlength' bl = qrkd_hat' :- halfbandlength' revb = J(1,length(qlist),1) revqlist = J(1,length(qlist),1) for( idx = 1 ; idx <= length(qlist); idx++ ){ revb[idx] = b[length(qlist)+1-idx] revqlist[idx] = qlist[length(qlist)+1-idx] } bbb95 = b[1] \ bu' \ revb' \ bl' \ b[length(qlist)] qqq95 = qlist[1] \ qlist' \ revqlist' \ qlist' \ qlist[length(qlist)] st_matrix(bname, b) st_matrix(Vname, V) st_numscalar(nname, n) st_numscalar(hname, h) st_matrix(cbname, (bbb95,qqq95)) st_matrix(qname, qlist) st_matrix(blname, bl) st_matrix(buname, bu) //////////////////////////////////////////////////////////////////////////// // Console output printf("\n{hline 78}\n") printf( "Number of observations: n = %f\n", n) printf( "The kink location of the running variable: k = %f\n", cut) printf( "The derivative of policy function on the left: b'(k-) = %f\n", bpleft) printf( "The derivative of policy function on the right: b'(k+) = %f\n", bpright) printf( "{hline 78}\n") printf( "Quantile QRKD [%2.0f%% Unif. Conf. Band]\n",100*cover) printf( "{hline 48}\n") for( idx = 1 ; idx <= length(qlist) ; idx++ ){ printf(" %5.3f %10.5f %10.5f %10.5f\n",qlist[idx],b[idx],bl[idx],bu[idx]) } printf( "{hline 78}\n") printf( "Test of the hypothesis that QRKD=0 for all quantiles: p-value = %4.3f\n", pvalueWS) printf( "Test of the hypothesis that QRKD is constant across quantiles: p-value = %4.3f\n", pvalueWH) printf( "{hline 78}\n") } end ////////////////////////////////////////////////////////////////////////////////