////////////////////////////////////////////////////////////////////////////////
// STATA FOR  Kato, K. & Sasaki, Y. (2018): Uniform Confidence Bands in
//            Deconvolution with Unknown Error Distribution. Journal of
//            Econometrics 207 (1), 129-161.
//
// Use it when you consider repeated measurements x1 and x2 for the true latent
// variable x with measurement errors e1 = x1 - x and e2 = x2 - x. The command
// draws a uniform confidence band for the density function fx of x.
////////////////////////////////////////////////////////////////////////////////
program define dkdensity, rclass
    version 14.2

    syntax varlist(min=2 numeric) [if] [in] [, numx(real 20), domain(real 2) cover(real 0.95), tp(real 0.2)]
    marksample touse

    gettoken depvar indepvars : varlist
    _fv_check_depvar `depvar'
    fvexpand `indepvars'
    local cnames `r(varlist)'

    tempname N CB x fx CBl CBu

		mata: estimate("`depvar'", "`cnames'", `numx', `domain', ///
					   `tp', `cover', "`touse'", "`N'", "`CB'", ///
					   "`x'", "`fx'", "`CBl'", "`CBu'")

	_matplot `CB', connect(1) noname ytitle("f{sub:X}(x)") xtitle("x") recast(line) name(UCB, replace)

	return scalar N = `N'
	return matrix upper = `CBu'
	return matrix lower = `CBl'
	return matrix fx = `fx'
	return matrix x = `x'
    return local cmd "dkdensity"
	di "*dkdensity command is based on Kato, K. & Sasaki, Y. (2018): Uniform Confidence"
	di "Bands in Deconvolution with Unknown Error Distribution. Journal of Econometrics"
	di "207 (1), 129-161."
end





mata:
////////////////////////////////////////////////////////////////////////////////
// Flat Top Kernel (Politis and Romano, 1999)
void phi_k(u, kout){
	k_b = 1.00
	k_c = 0.05
	abs_u = abs(u)
	kout = exp( -1 :* k_b :* exp(-1 :* k_b :/ ((abs_u :- k_c):^2 + 0.000001)) :/ ((abs_u :- 1):^2 + 0.000001) )
	kout = 1 :* (abs_u :<= k_c) :+ kout :* (k_c :< abs_u)
	kout = kout :* (abs_u :< 1) :+ 0 :* (1 :<= abs_u)
}
////////////////////////////////////////////////////////////////////////////////
// Estimation of QRKD
void estimate( string scalar x1v,     string scalar x2v,
			   real scalar numx,      real scalar domain,
			   real scalar tuning,	  real scalar cover,
			   string scalar touse,   string scalar nname,
			   string scalar cbname,
			   string scalar xname,   string scalar fxname,
			   string scalar cblname, string scalar cbuname)
{
    x1      = st_data(., x1v, touse)
    x2      = st_data(., x2v, touse)
    n      = length(x1)

	// Normalize the Location of x1 to zero
	mean_x1 = mean(x1)
	sd_x1 = ( mean( x1 :* x2 ) :- mean(x1):*mean(x2) )^0.5 // variance(x1)^0.5
	x1 = ( x1:-mean_x1 ) :/ sd_x1
	x2 = ( x2:-mean_x1 ) :/ sd_x1

	// Define y and eta
	y = ( x1 :+ x2 ) :/ 2
	eta = ( x1 :- x2 ) :/ 2

	// xlist
	xlist = ( ( (1..numx) :- ( (1+numx)/2 ) ) :/ ( (numx-1)/2 ) ) :* domain

	// tlist
	tlist = (-20..20) / 10
	tlist_interval = tlist[2] - tlist[1]

	// Bandwidth
	h = tuning

	////////////////////////////////////////////////////////////////////////////
	// Compute the Fourie Transforms
	phi_K = tlist :* (0+0i)
	for( idx = 1 ; idx <= length(tlist) ; idx++ ){
		t = tlist[idx]
		phi_k(t, phi_k_val)
		phi_K[idx] = phi_k_val
	}

	phi_e = tlist :* (0+0i)
	for( idx = 1 ; idx <= length(tlist) ; idx++ ){
		t = tlist[idx]
		phi_e[idx] = sum( exp( (0+1i) :* t :/ h :* eta ) ) / n
	}
	//phi_K',phi_e'

	////////////////////////////////////////////////////////////////////////////
	// Compute Density and Variance
	fx = xlist :* 0
	sigma2 = xlist :* 0
	K_hat = J(n,length(xlist),0+0i)

	for( idx = 1 ; idx <= length(xlist) ; idx++ ){
		x = xlist[idx]

		for( jdx = 1 ; jdx <= n ; jdx++ ){
			K_hat[jdx,idx] = sum( exp( (0-1i) :* tlist :* (x :- y[jdx]):/h ) :* phi_K :/ phi_e ) :* tlist_interval :/ (2 * 3.14159265359)
		}
		sum_k_hat = sum( K_hat[.,idx] )
		sum_k_hat_sq = sum( K_hat[.,idx]:^2 )

		fx[idx] = Re( sum_k_hat / (n*h) )
		sigma2[idx] = Re( sum_k_hat_sq / n - ( sum_k_hat / n )^2 )
	}
	//fx' :/ sd_x1, sigma2' :/ sd_x1:^2

	////////////////////////////////////////////////////////////////////////////
	// Multiplier Bootstrap
	num_boot = 4000
	supnormZxi = J(1,num_boot,0)

	for( jdx = 1 ; jdx <= num_boot ; jdx++ ){
		xi = invnormal(uniform(n,1))
		Zxi = xlist * 0
		for( idx = 1 ; idx <= length(xlist) ; idx++ ){
			Zxi[idx] = Re( sum( xi :* ( K_hat[.,idx] :- sum( K_hat[.,idx] ) :/ n ) ) / ( sigma2[idx]:^0.5 * n^0.5 ) )
		}
		supnormZxi[jdx] = max(abs(Zxi))
	}
	critical_value = supnormZxi[trunc(cover*num_boot)]
	//critical_value

	////////////////////////////////////////////////////////////////////////////
	// Confidence Bands
	CBl = fx :- sigma2:^0.5 :/ ( n^0.5 * h ) :* critical_value
	CBu = fx :+ sigma2:^0.5 :/ ( n^0.5 * h ) :* critical_value
	//CBl', fx', CBu'

	////////////////////////////////////////////////////////////////////////////
	// Final Output
	final_xlist = xlist' :* sd_x1 :+ mean_x1
	final_fx = fx' :/ sd_x1
	final_CBl = CBl' :/ sd_x1
	final_CBu = CBu' :/ sd_x1

	//final_xlist, final_fx, final_CBl, final_CBu

	////////////////////////////////////////////////////////////////////////////
	// Display Results
	for( idx = 1 ; idx <= length(xlist) ; idx++ ){
	if( trunc((idx-1)/10) == (idx-1)/10 ){
		printf("{hline 64}\n")
		printf("         x           f(x)     [%2.0f%% Uniform Confidence Band]\n", 100*cover)
		printf("{hline 64}\n")
	}
	printf("%10.3f     %10.3f     %10.3f     %10.3f\n", final_xlist[idx], final_fx[idx], final_CBl[idx], final_CBu[idx])
	}
	printf("{hline 64}\n")

	////////////////////////////////////////////////////////////////////////////
	// Set
	rev_final_fx = J(length(xlist),1,1)
	rev_final_xlist = J(length(xlist),1,1)
	for( idx = 1 ; idx <= length(xlist); idx++ ){
	    rev_final_fx[idx] = final_fx[length(xlist)+1-idx]
	    rev_final_xlist[idx] = final_xlist[length(xlist)+1-idx]
	}
	BBB = final_fx[1] \ final_CBu \ rev_final_fx \ final_CBl \ final_fx[length(xlist)]
	XXX = final_xlist[1] \ final_xlist \ rev_final_xlist \ final_xlist \ final_xlist[length(xlist)]
    st_matrix(cbname, (BBB,XXX))

    st_numscalar(nname, n)
	st_matrix(xname, final_xlist)
	st_matrix(fxname, final_fx)
	st_matrix(cblname, final_CBl)
	st_matrix(cbuname, final_CBu)
}
end
////////////////////////////////////////////////////////////////////////////////