@ Computation of p-value for the Andrews (1993) stability test given the critical value C. The algorithm is based on    @
@ Anatolyev and Kosenok (2011) "Another numerical method of finding critical values for the Andrews stability test"     @ 

//************ SPECIFY YOUR PARAMETERS HERE ************
P1 = 0.05;      @ left end of boundary         		@
P2 = 0.95;      @ right end of boundary         	@
Ka = 3;         @ degrees of freedom for Bessel bridge 	@
C = 15.53; 	@ critical value of boundary          	@
//******************************************************

CalcAlpha( P1, P2, C, Ka ); end;

@ Calculates the size for the boundary c*r*(1-r) on the segment [P1,P2]   @
proc (1) = CalcAlpha( P1, P2, c, k );
local m, eps, f, hi, gr, bnd, ff, fb, hw, w, bw, v0, wm, bwm, v1, b, b1, s, h, Kii0, Kii1, Kii, x0, x1, y, ds, a;

m = 501;
eps = 0.000001;
f = zeros(m,1);
hi = 0.001;
gr = seqa(0,1/m,m);
gr = gr*(P2-P1-hi)/gr[m]+P1+hi;
gr = P1|gr;
bnd = c*gr.*(1-gr);

ff = zeros(5,5);
fb = zeros(5,1);
for i (1,5,1);
hw = i*hi/5;
w = gr[1]+hw; bw = c*w*(1-w);
v0 = Prs(gr[1],w,bnd[1],bw,k);
wm = gr[1]+hw/2; bwm = c*wm*(1-wm); // boundary
v1 = Prs(wm,w,bwm,bw,k);
b = -v0*hw^(1/2)+hw^(1/2)*2^(1/2)*v1;
b1 = (2^(1/2)*v0-v1)*2^(1/2)/hw^(1/2);
fb[i] = VerI(gr[1],w,bnd[1],bw,k); 

ff[i,1] = b*Pi+b1*Pi*hw/2;
ff[i,2] = 2*b*hw^(1/2)+(2/3)*b1*hw^(3/2);
ff[i,3] = (1/8)*hw*Pi*(4*b+b1*hw);
ff[i,4] = (4/3)*b*hw^(3/2)+(4/15)*b1*hw^(5/2);
ff[i,5] = (1/16)*hw^2*Pi*(6*b+b1*hw);
endfor;

a = inv(ff)*fb;
s = sqrt(hi);
f[1] = a[1]/s+a[2]+a[3]*s + a[4]*s^2+ a[5]*s^3;
h = gr[3]-gr[2];

for i (2,m,1);
  f[i] = VerI(gr[1],gr[i+1],bnd[1],bnd[i+1],k);
  w = gr[i+1] - 2*eps; bw = c*w*(1-w);
  Kii0 = Prs(w,gr[i+1],bw,bnd[i+1],k)*sqrt(2*eps);
  w = gr[i+1] - eps; bw = c*w*(1-w);
  Kii1 = Prs(w,gr[i+1],bw,bnd[i+1],k)*sqrt(eps);
  Kii = 2*Kii1-Kii0;

  v0  = Prs(gr[1],gr[i+1],bnd[1],bnd[i+1],k); // integral on the first segment
  v1  = Prs(gr[2],gr[i+1],bnd[2],bnd[i+1],k);
  
  fb[1] = (4/3)*hi^(1/2)*v0+(2/3)*hi^(1/2)*v1;
  fb[2] = (1/2)*v0*hi+(1/2)*hi*v1;
  fb[3] = (2/5)*v1*hi^(3/2)+(4/15)*v0*hi^(3/2);
  fb[4] = (1/3)*hi^2*v1+(1/6)*v0*hi^2;
  fb[5] = (2/7)*hi^(5/2)*v1+(4/350)*v0*hi^(5/2);

  s = a'*fb; 
  
  v1 = Prs(gr[2],gr[i+1],bnd[2],bnd[i+1],k)*f[1]*sqrt(gr[i+1]-gr[2]);
  for j (1,i-2,1);
    v0 = v1;
    v1 = Prs(gr[j+2],gr[i+1],bnd[j+2],bnd[i+1],k)*f[j+1]*sqrt(gr[i+1]-gr[j+2]);
    x0 = gr[j+1]; x1 = gr[j+2]; y = gr[i+1];
    ds = -(2/3)*(-3*v0*x1*(y-x0)^(1/2)-2*v1*y*(y-x0)^(1/2)+2*v1*x0*(y-x0)^(1/2)+2*v0*y*(y-x0)^(1/2)+v0*x0*(y-x0)^(1/2)
         +2*v0*(y-x1)^(1/2)*x1+2*v1*y*(y-x1)^(1/2)+v1*(y-x1)^(1/2)*x1-3*v1*x0*(y-x1)^(1/2)-2*v0*y*(y-x1)^(1/2))/(x1-x0);
    s = s + ds;
  endfor;
    v0 = Prs(gr[i],gr[i+1],bnd[i],bnd[i+1],k)*f[i-1]*sqrt(h); 
    v1 = Kii;
    x0 = gr[i]; x1 = gr[i+1];
    f[i] = -(1/4)*(2*(x1-x0)^(1/2)*v0-3*(f[i]-s))/(x1-x0)^(1/2)/Kii;
  
endfor;
s = sqrt(hi);
retp( sumc(f)*h-0.5*f[1]*h-0.5*f[m]*h + 2*s*a[1] + a[2]*s^2 + a[3]*(2/3)*s^3 + a[4]*s^4/2 + a[5]*(2/5)*s^5 + cdfchic(c,k) );
endp;


@ Calculates the RHS of the integral                            @
proc (1) = VerI( s, r, x, y, k );
local xx, v0, v1, t, f, NIC;
NIC=100;
xx = x;
v1 = Pr(s,xx,k)*Prs(s,r,xx,y,k);
xx = x/2;
v0 = Pr(s,xx,k)*Prs(s,r,xx,y,k);
do while v0 < v1*1.e-16;
xx = 0.5*(x+xx);
v0 = Pr(s,xx,k)*Prs(s,r,xx,y,k);
endo;
xx = xx-(x-xx);
t = (x-xx)*((cos(Pi*(seqa(0,1,NIC)+0.5)/NIC)+1)/2)+xx;
f = t;
for i (1,NIC,1);
f[i] = Pr(s,t[i],k)*Prs(s,r,t[i],y,k); 
endfor;
retp( (x-xx)*CCQ100(f) );
endp;

@ Calculates the unconditional pdf of the integral equation                @
proc (1) = Pr( r, y, k );
retp( (2*r*(1-r))^(-k/2)*y^(k/2-1)*exp(-y/(2*r*(1-r)))/GAMMA(k/2) );
endp;

@ Calculates the conditional pdf under the integral equation                @
proc (1) = Prs( s, r, x, y, k );
local wsr;
wsr = (1-s)/((1-r)*(r-s));
retp( wsr*DCHINC(wsr*y, k, sqrt(x/wsr)/(r-s) ) );
endp;

@ Calculates density of chi-square noncentral distribution                  @
proc (1) = DCHINC( z, v, d );
local x1, x2, CMUL, N, M, f, r, t;
N = 100;  M=49;
if v < 2; 
retp( (exp(-0.5*(d-sqrt(z))^2)+exp(-0.5*(d+sqrt(z))^2))/(2*sqrt(2*Pi*z)) );
else;
CMUL = 2^((v-3)/2)*gamma((v-1)/2)*sqrt(2*Pi);
r = minc(Pi|(8*sqrt((v)/(sqrt(z)*d))) );
t = r*(cos(Pi*(seqa(0,1,N)+0.5)/N)+1)/2;
x1 = d-sqrt(z)*cos(t);
x2 = sqrt(z)*sin(t);
f = (exp(-0.5*(x1^2+x2^2)).*x2^(v-2)/CMUL)/2;
retp( r*CCQ100(f) );
endif;
endp;

@ Clenshaw-Curtis Quadrature with 100 points    @
proc (1) = CCQ100( f );
local Tj;
let Tj = {
0.000215377999909553 
0.000757443464203025 
 0.00122102983079868 
 0.00173191300275900 
 0.00220701989824522 
 0.00270557138012544 
 0.00318144060915189 
 0.00367012706637387 
 0.00414234327043413 
 0.00462083278245196 
 0.00508645475626408 
 0.00555362672387725 
 0.00601024529199005 
 0.00646469739786849 
 0.00691015877417317 
 0.00735038551432303 
 0.00778269010988163 
 0.00820716102926611 
 0.00862442215655300 
 0.00903162225126352 
 0.00943204903015457 
 0.00982050257544876 
  0.0102023938215584 
  0.0105706800463154 
  0.0109324235527214 
  0.0112791878922866 
  0.0116192624856046 
  0.0119432252171702 
  0.0122602042526497 
  0.0125601674479771 
  0.0128527230115735 
  0.0131275763197367 
  0.0133944837089402 
  0.0136432092639678 
  0.0138833514815563 
  0.0141050281116364 
  0.0143174001980853 
  0.0145112070459907 
  0.0146949201306023 
  0.0148601397555587 
  0.0150144247403034 
  0.0151504457474962 
  0.0152746565599429 
  0.0153809757887377 
  0.0154745921560889 
  0.0155508164481967 
  0.0156134461560861 
  0.0156592937182519 
  0.0156906743271761 
  0.0157059756982723 
  0.0157059756982723 
  0.0156906743271761 
  0.0156592937182519 
  0.0156134461560861 
  0.0155508164481967 
  0.0154745921560889 
  0.0153809757887377 
  0.0152746565599429 
  0.0151504457474962 
  0.0150144247403034 
  0.0148601397555587 
  0.0146949201306023 
  0.0145112070459907 
  0.0143174001980853 
  0.0141050281116364 
  0.0138833514815563 
  0.0136432092639678 
  0.0133944837089402 
  0.0131275763197367 
  0.0128527230115735 
  0.0125601674479771 
  0.0122602042526497 
  0.0119432252171702 
  0.0116192624856046 
  0.0112791878922866 
  0.0109324235527214 
  0.0105706800463154 
  0.0102023938215584 
 0.00982050257544877 
 0.00943204903015457 
 0.00903162225126352 
 0.00862442215655300 
 0.00820716102926612 
 0.00778269010988162 
 0.00735038551432303 
 0.00691015877417318 
 0.00646469739786850 
 0.00601024529199006 
 0.00555362672387725 
 0.00508645475626408 
 0.00462083278245196 
 0.00414234327043413 
 0.00367012706637388 
 0.00318144060915189 
 0.00270557138012544 
 0.00220701989824522 
 0.00173191300275900 
 0.00122102983079868 
0.000757443464203027 
0.000215377999909556 
};
retp(Tj*f);
endp;