Re: st: sign test output

[Maarten Buis <maartenlbuis@gmail.com>]

On Fri, Jan 18, 2013 at 10:48 PM, Dirk Enzmann wrote:
> In this context, the following working paper by Mantalos might be
> interesting:
>
> http://hj.se/download/18.3bf8114412e804c78638000150/1299244445855/WP2010-8.pdf
>
> It should be possible to implement his JBCV(k1,k2) procedure in Stata and it
> would be interesting to see the results of including this test in the
> simulation.

What they do is in the bootstrap literature called an Achieved
Significance Level (ASL), there is an example in the manual entry for
-bootstrap- and I have given an example of that recently in this
thread (for a one sample t-test). Below I have written a small program
to compute the ASL for the Jarque-Bera test.

ASLs are computed with aid of a Monte Carlo experiment and thus
uncertain. This program also quantifies this uncertainty by displaying
a Monte Carlo confidence interval for the ASL, so if you were to call
-jb_asl- a 100 times you would expect to find an ASL within this
interval in 95 of these calls. If you are unhappy with the width of
this interval, than you can decrease it by specifying more
replications in the -reps()- option. Typically, I find 20,000 an
number that works well for this type of procedure; if the "true"
p-value is 0.05 than the ASL will on average be computed based on a
1000 replications where the simulated J-B statistic is larger than the
observed J-B statistic.

As to how it will perform in simulations I can tell two things: 1) it
will be slow as it is in essence a simulation run within a simulation,
and 2) the p-values will at best follow a discrete uniform
distribution rather than a continuous uniform distribution. I don't
think that point 2 is a big issue, as you can arbitrarily approximate
the continuous distribution by increasing the number of replications,
but for a simulation you would typically need to specify only a few
replications in order to deal with problem 1). So, when evaluating how
well this procedure works with a simulation one needs to take problem
2) into account.

*----------------------------------- begin program ----------------------------
*! version 1.0.0 20Jan2013 MLB
program define jb_asl
syntax varname [if] [in],             ///
    [                                 ///
	reps(numlist max=1 integer >0)    ///
	nodots                            ///
	mcci(cilevel)                     ///
	]
	
	
	if "`reps'" == "" {
		local reps = 1000
	}
	tempname scalar jb jbp m sd asl jb_sim count lb ub alph
	tempvar x
	marksample touse
	qui sum `varlist' if `touse', detail
	scalar `jb'  = (r(N)/6) * ///
				   (r(skewness)^2 + 1/4*(r(kurtosis) - 3)^2)
	scalar `jbp' = chi2tail(2,`jb')
	scalar `m'   = r(mean)
	scalar `sd'  = r(sd)
	qui gen `x' = .
	scalar `asl' = 0
	scalar `count' = 0
	if "`dots'" == "" {
		_dots 0, title(computing ASL) reps(`reps')
	}
	forvalues i = 1/`reps' {
		capture {
			replace `x' = rnormal(`m',`sd') if `touse'
			sum `x' if `touse', detail
			scalar `jb_sim' = (r(N)/6) * ///
					   (r(skewness)^2 + 1/4*(r(kurtosis) - 3)^2)
			scalar `asl' = `asl' + (`jb_sim' > `jb')
		}
		if !_rc {
			scalar `count' = `count' + 1
		}
		if "`dots'" == "" {
			_dots `i' `=_rc>0'
		}
	}
	scalar `alph' = (100-`mcci')/200
    scalar `lb' = invbinomial(`count',`asl',   `alph')
    scalar `ub' = invbinomial(`count',`asl', 1-`alph')
	scalar `asl' = `asl'/`count'
	di as txt "Jarque-Bera statistic:               " as result %6.2f `jb'
	di as txt "asymptotic p-value:                  " as result %6.4f `jbp'
	di as txt "achieved significance level (ASL):   " as result %6.4f `asl'
	di as txt "`mcci'% Monte Carlo CI for ASL: {col 38}[" _c
	di as result %6.4f `lb' as txt ", " as result %6.4f `ub' _c as txt "]"
end
*---------------------------------- end program --------------------------------

-- Maarten

---------------------------------
Maarten L. Buis
WZB
Reichpietschufer 50
10785 Berlin
Germany

http://www.maartenbuis.nl
---------------------------------
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/