McCrary density diagnostics - testing for heaping at round numbers - STATA Commands

[Ammar Hussain]

Hi,

My running variable is the price ratio, a continuous variable that is equal to the average quarter share price / recent offering price. In my case, the cutoff point is "1" (when the average quarter share price becomes equal to the recent offering price).


For testing the manipulation in my running variable, I used the following McCrary test command in STATA:

rddensity PriceRatio, plot c(1) p(1)

But, my supervisor suggested: The denominator (last offer price) is endogenous to earlier timing and may mechanically create mass around ratios near 1. So, also testing for heaping at round numbers, which enhances the quality of the analysis.


MY QUESTION(s): 1) What is testing heaping at round numbers, and why do we test this?
2) Please suggest to me the actual codes/commands that I used to test the heaping at the round number based on the above case.



Thanks

[George Ford]

The analysis of heaping is often visual (spikes in histogram) or by tabulations.

https://jmparman.people.wm.edu/stata...-examples.html

This is a proposed test:

https://www.tandfonline.com/doi/pdf/...64760120074960

[Ammar Hussain]

George Ford , could you please share the STATA commands for testing heaping at round numbers through the rddensity test for the above-mentioned running variable and cutoff value? I searched and read much literature, but did not understand for my case.

Thanks in Advance for your valuable time and precious guidance.

[George Ford]

I don't have the code, but it shouldn't be hard to code that test.

[Ammar Hussain]

Can anyone else kindly respond to it?

[George Ford]

If you want someone to do your work for you, then you'll probably need to pay them.

[George Ford]

This is a rewrite of nipnTK::ageHeaping() from R.

y has no heaping
yh has heaping.

Both are tested.

clear all

set obs 1000

g y = runiformint(10,20)
g yh = y
replace yh = 15 if runiform()<0.1
hist yh
unique yh

capture program drop testheaping
program define testheaping, rclass
version 15.1
// Usage: ageheaping agevar, [divisor(#)]
syntax varname(numeric) [if] [in], [DIVISOR(integer 12)]

if (`divisor' <= 1) {
di as err "divisor() must be an integer > 1"
exit 198
}

marksample touse, strok

// Remainder
tempvar rem
quietly gen double `rem' = mod(`varlist', `divisor') if `touse' & !missing(`varlist')

// Total N (nonmissing in sample)
quietly count if !missing(`rem')
local N = r(N)
if (`N' == 0) {
di as err "No non-missing observations in sample."
exit 2000
}

// Build a full 0..divisor-1 remainder table (including zero-count cells)
tempname counts pct
matrix `counts' = J(`divisor', 1, 0)

forvalues k = 0/`= `divisor' - 1' {
quietly count if `rem' == `k'
matrix `counts'[`= `k' + 1', 1] = r(N)
}

// Percent table (one decimal, like R code's round(..., 1))
matrix `pct' = (`counts' / `N') * 100
forvalues i = 1/`divisor' {
matrix `pct'[`i',1] = round(`pct'[`i',1], .1)
}

// Chi-square goodness-of-fit to uniform over remainders
// expected count in each cell:
scalar exp = `N' / `divisor'
scalar chi2 = 0
forvalues i = 1/`divisor' {
scalar obs = `counts'[`i',1]
// include zero-count categories (as in nipnTK fullTable)
scalar chi2 = chi2 + ((obs - exp)^2) / exp
}
scalar df = `divisor' - 1
scalar p = chi2tail(df, chi2)

// Name rows nicely: remainder 0..divisor-1
local rnames
forvalues k = 0/`= `divisor' - 1' {
local rnames `"`rnames' `k'"'
}
matrix rownames `counts' = `rnames'
matrix rownames `pct' = `rnames'
matrix colnames `counts' = count
matrix colnames `pct' = pct

// Display similar to R print()
di _n as txt "{tab}Heaping Analysis" _n
di as txt "data:{tab}Remainder of `varlist' / `divisor'"
di as txt "X-squared = " %10.4f chi2 ///
as txt ", df = " %6.0f df ///
as txt ", p-value = " %9.6f p
di _n as txt "Remainder distribution (counts):"
matlist `counts', names(rows) format(%12.0f)
di _n as txt "Remainder distribution (percent):"
matlist `pct', names(rows) format(%9.1f)

// Return results
return scalar N = `N'
return scalar X2 = chi2
return scalar df = df
return scalar p = p
return matrix tab = `counts'
return matrix pct = `pct'
end

testheaping y, divisor(11)

testheaping yh, divisor(11)