/*
** TESTING MEAN BY BOOTSTRAP **
**    X IID FROM CHI^2(5)    **
**   MONTE CARLO EXERCISE    **
**    PREPARED BY S.C. AHN   */


@ CLEAN UP MEMORY @

    new ;

@ OPEN OUTPUT FILE @

    output file = bmeanmon.out reset ;

@ SET-UP: CAN CHANGE @

    sd    = 1    ; @ initial seed number for random number generation @
    tt    = 100  ; @ number of observation @
    nboot = 1000 ; @ number of bootstrap samples @
	iter  = 1000 ; @ number of simulations @
	mu    = 5    ; @ value of mu under H_0 @

@ INITIAL MATRICES: DO NOT CHANGE @

	abias = zeros(iter,1) ;
	bbias = zeros(iter,1) ;
	amse  = zeros(iter,1) ;
	bmse  = zeros(iter,1) ;
	arej  = zeros(iter,3) ;
	brej  = zeros(iter,3) ;

/* STARTING MONTE CARLO */

j = 1 ; do while j <= iter ;


@ X IID FROM CHI(KK): CAN CHANGE @

	sd = j + 100                  ; @ updating seed number @                      
    kk = 5                        ; @ degrees of freedom @
    xx = sumc(rndns(kk,tt,sd)^2)  ; @ tt*1 data vector @

/* DO NOT CHANGE FROM HERE */

@ COMPUTING SAMPLE MEAN AND SAMPLE VARIANCE @

    xb  = meanc(xx)   ; @ sample mean     @
    xs2 = stdc(xx)^2  ; @ sample variance @

	abias[j] = xb - kk    ;
	amse[j] = (xb - kk)^2 ;

@ ASYMPTOTIC CHI_SQUARE WALD TEST FOR H_0: MU = TRUE MU @

    awt  = tt*(xb-mu)^2/xs2        ; @ chi^2 Wald stat.  @

	if awt > 6.63490 ; test1 = 1 ; else ; test1 = 0 ; endif ;
	if awt > 3.84146 ; test5 = 1 ; else ; test5 = 0 ; endif ;
	if awt > 2.70554 ; test0 = 1 ; else ; test0 = 0 ; endif ;

	arej[j,.] = test1~test5~test0 ;	

/* BOOTSTRAP WITH PROB = 1/T FOR ALL OBSERVATIONS */

bootxb = zeros(nboot,1) ;
bootwt = zeros(nboot,1) ;

i = 1 ; do while i <= nboot ;

    cboot = ceil( tt*rndus(tt,1,sd) ) ; @ choosing bootstrap data points @
    bsam  = xx[cboot,.] ;

    bxb   = meanc(bsam)  ; @ boot-sample mean     @
    bxs2  = stdc(bsam)^2 ; @ boot-sample variance @

    bwt       = tt*(bxb-xb)^2/bxs2 ; @ boot-chi^2 Wald stat. @
    bootxb[i] = bxb                ; @ saving boot-means     @
    bootwt[i] = bwt                ; @ saving boot-stat      @

i = i + 1 ; endo ;

@ BIAS @

    bias = meanc(bootxb) - xb ;

@ BIAS-CORRECTED ESTIMATE @

    bcxb = xb - bias ;

	bbias[j] = (bcxb - kk)   ;
	bmse[j]  = (bcxb - kk)^2 ;

@ BOOTSTRAP CRITICAL VALUES @

    bootwt = sortc(bootwt,1) ;
    tc     = nboot           ;
    bootc  = bootwt[ceil(tc*0.99)]|bootwt[ceil(tc*0.95)]|bootwt[ceil(tc*0.9)];

	if awt > bootc[1] ; test1 = 1 ; else ; test1 = 0 ; endif ;
	if awt > bootc[2] ; test5 = 1 ; else ; test5 = 0 ; endif ;
	if awt > bootc[3] ; test0 = 1 ; else ; test0 = 0 ; endif ;

	brej[j,.] = test1~test5~test0 ;

j = j + 1 ; endo ;


/* REPORTING RESULTS */

format /rd 12,4 ;

"BIASES IN POINT ESTIMATES " ;
"" ;
"   SAMPLE MEAN  BIAS-COR. MEAN" ;
meanc(abias) meanc(bbias) ;
"";
"REJECTION RATE OF ASYMPTOTIC WALD TEST FOR H_0: MU = " mu ;
"";
"      at 1%        at 5%        at 10%";
meanc(arej)' ;
"" ;
"BOOTSTRAP WALD TEST FOR H_0: MU =  " mu ;
"";
"      at 1%        at 5%        at 10% ";
meanc(brej)' ;

output off ;