/**********        Bradley & McClelland Misspecification Test   **********************

This is the Gauss program for the test described in the paper "An
Improved Nonparametric Test for Misspecification of Functional Form" by
Ralph Bradley and Robert McClelland.  It is strongly recommended that
users read that paper before using this test. A postscript version of
the paper is available at
http://stats.bls.gov/cgi-bin/OREPrintAbstract/econ.1993.bradley_ralph

Out test is designed to detect specification errors in the first
conditional moment of the regression (i.e. the functional form of the
regression) on iid data.  In our paper We show that it is consistent
against all alternatives from the null hypothesis of correct
specification (excepting a parameter set of measure zero).  We also
show that it is the asymptotically most powerful test in the sense
that the plim of the test maximizes plim W/n, where W is the set of
conditional moment tests with chi-square 1 distributions and n is the
sample size.

The basic intuition is that we wish to choose a function that maximizes
u_hat*g(x)/V, where u_hat is the estimated residuals from a nonlinear
least squares regression, g(x) is any x-measurable function and V is
a consistent estimate of the variance, when the null hypothesis is false
and equal to zero when true.  By the Rao-Blackwell theorem, the function
g(x) that maximizes this function is E(u_hat|x).  If the user has
parametric information about the nature of this function, it can be used in
the test.  If the researcher does not, as is often the case, have any
information, a nonparametric estimator should be used.  We implement
this by using a kernel density estimate for the nonparametric estimator.

This test can be implemented two ways.  If the user has parametric
information about E(u_hat|x), then use test(u_hat,g_hat,x), where
g_hat is the estimate of E(u_hat|x) and x is the vector of partial
derivitives of the dependent variable with respect to the regressors.
Each variable in the vector must be continuous.  If the model is
linear, then x is obviously just the set of regressors themselves.

If the user has no information then call ktest(x, u_hat, Nobs2, iseed),
where Nobs2 is the number of observations for the bootstrapped sample
(necessary to reduce bias in the kernel density estimate) and iseed is
a seed number for the random-number generator for the bootstrapped
sample.  Nobs2 should be O(n), but there is otherwise no restrictions
put on it.

As mentioned above, the result is a random variable that should have a
chi-squared one distribution under the null.  Also, be warned that this
test is very slow, due to its thoroughness.

Everything in this program is in the public domain, but if you use the
test please cite us and the paper as the source.  Also, we bear no
responsibility for any problems associated with the use of this program.

 Rob McClelland      McClelland_R@dcgate.bls.gov     9/15/94
***************************************************************************/



@***************************************************************************
            KERNEL ESTIMATION SECTION
    This section lists the procedures for kernel estimation of a
    regression function
***************************************************************************@

proc nahm(z) ;
  @ one possible kernel @
  local a,b,c,d ;

  d = sumc((-z.*z/2)') ;
  c =  abs(d) ;
  d = d.*(c .lt 35) ;
  b = (c .lt 32).*exp(d) ;
retp(b) ;
endp ;



proc  kerreg(y,x,n,k,w,xb,&f1) ;
/* kernel regression. Does not correct for bias */
/* Guts of the instrumental kernal   y,x  = sample */
/*                                   yinst = instrumental vector */
/*                                   n = sample size       */
/*                                   k = dim x             */
/*                                   w = window width      */
/*                                   xb = value to compute regression */
/*                                   &f1 = kernel function */
/*                                   xd = alternative instuments */


 local f1:proc ;
 local a,b,c,ck,xs,xf,kern, i2,i3,xmean,xstd,xvar,yinm,yins;
 local sumnum,sumdem,ye1r,s1mdem,s1mnum,s2mdem,s2mnum,f;

 sumnum = 0.0 ;
 sumdem = 0.0 ;

 xs = x - xb;

 /* print  x;"xs";  xs ;"xb";  xb ;  */

 @ Normalization making c identical up to a linear transformation @
 @ See Bierens p.125 eq 5.1.0                                     @

 xmean = meanc(xs)^2;
 xstd  = sqrt(stdc(xs).^2*(n/(n-1))^2 + xmean);

 c = (xs./xstd')./w;

 flow:

 ck = f1(c) ;


 sumnum = y'*ck ;
 sumdem = sumc(ck) ;

 if (sumdem < 0.00001) ;
    sumdem = 0.0 ;
    c = inv(1.25).*c ;
    goto flow ;
 endif;

 kern = sumnum / sumdem ;

 retp(kern) ;
 endp ;



proc (1) = ukerreg(y,x,c,xb,&f1,delta) ;
  /* Jacknife Kernel to purge the bias.  Note -- c is not window.
 The window is wv=cn^(1/(1+k)) */

local f1:proc ;
local y1,y2,y3,wv,adj,n,k ;
k=cols(x);
n=rows(x);

wv = c*inv(n^(inv(4+k)));

y1 = kerreg(y,x,n,k,wv,xb,&f1) ;

wv = c*inv(n^(delta/(4+k)));
y2 = kerreg(y,x,n,k,wv,xb,&f1) ;

adj= inv(n^(2*(1-delta)/(k+4)));
y3 = (y1 - adj*y2)/(1-adj) ;

 retp(y3~y1);
 endp  ;


@***********************************************************************
            END OF KERNEL ESTIMATION SECTION
************************************************************************@

@***************************************************************************
            OPTIMIZATION SECTION
    This section lists the procedures for optimization using the golden
    rule routine and cross-validation minimizing out-of-sample MSE
***************************************************************************@


proc cross(y,x,n,k,cc) ;
@ selects optimal window width by minimizing out-of-sample MSE @
@ cc is the window width @

local er,n1,ye  ;
local yn,xn,xb ;
local crit,d,xm,zx  ;
local be,by ;

er = ones(n,1)  ;
let ye[1,2] = 0;

d=.5;
xn = x[2:n,.];
yn = y[2:n];
xb = x[1,.] ;
ye =  ukerreg(yn,xn,cc,xb,&nahm,d) ;
er[1,1] = (y[1] - ye[1,1])^2 ;
xn = x[1:n-1,.];
yn = y[1:n-1];
xb = x[n,.] ;
ye =  ukerreg(yn,xn,cc,xb,&nahm,d) ;
er[n,1] = (y[n] - ye[1,1])^2 ;

n1=2;

do while n1 < n;

    xn = x[1:n1-1,.]|x[n1+1:n,.] ;
    yn = y[1:n1-1]|y[n1+1:n] ;
    xb = x[n1,.] ;

    ye =  ukerreg(yn,xn,cc,xb,&nahm,d) ;

    /* ye = kerreg(yn,xn,n-1,k,w,xb,&nahm) ; */
    er[n1,1] = (y[n1] - ye[1,1])^2 ;


    /*  print n1  "y" y[n1,1] "ye" ye zx[n1,1] "e2" er[n1,1] ; */
    n1 = n1 + 1 ;
endo;


crit = meanc(er) ;

retp(crit) ;
endp ;



proc(2) =  search(L_b,M_p,U_b,tol,y,x,n,k)  ;
/* Golden Rule Minimization */
local  ratio,window,MSE,MSE_1,MSE_2,c,w_1,w_2,w_3,w_4;

 ratio = 0.618;
 c = 1.0-ratio;
 w_1 = L_b;
 w_4 = U_b;

 w_3=M_p+(abs(U_b-M_p) .gt abs(M_p-L_b))*c*(U_b-M_p);
 w_2=M_p-(abs(U_b-M_p) .le abs(M_p-L_b))*c*(M_p-L_b);

 MSE_1 = cross(y,x,n,k,w_2) ;
 MSE_2 = cross(y,x,n,k,w_3);

  do WHILE abs(w_4-w_1) > tol*(abs(w_2)+abs(w_3));
  If MSE_2 < MSE_1  ;
     w_1 = w_2;
     w_2 = w_3;
     w_3 = ratio*w_2+c*w_4;
     MSE_1 = MSE_2;
     MSE_2 = cross(y,x,n,k,w_3) ;
  Else;
     w_4 = w_3;
     w_3 = w_2;
     w_2 = ratio*w_3+c*w_1;
     MSE_2 = MSE_1;
     MSE_1 = cross(y,x,n,k,w_3) ;
  Endif  ;
  Endo;

MSE= minc(MSE_1|MSE_2);
window = (MSE_1 .lt MSE_2).*w_2 + (MSE_1 .ge MSE_2).*w_3;

retp(MSE,window) ;
ENDp;



@************************************************************************
            END OF OPTIMIZATION SECTION
************************************************************************@


@*************************************************************************
            TEST SECTION
   This section lists the actual procedures called when invoking the test
   and the bootstrapping procedure as well.
**************************************************************************@

proc boot(x,uhat,Nobs2,iseed);

/* boot creates a bootstrapped sample and computes the optimal window
   width for that sample using a golden rule search routine.  It then
   estimates a jack-knived kernel regression for every point in the
   original data set, using a separate bootstrap for each point. To change
   the kernel from Gaussian to one the user supplies, change the name
   &nahm to the name of the supplied kernel.   The upper and lower
   bounds on the window width are necessary to guarantee that the statistic
   doesn't degenerate.  The bounds used here are just suggestions.  It
   is important that the bounds increase to eventually cover all positive
   real numbers as the number of observations increases to infinity. */

local bootidx, Nobs, x_boot, uhatboot, ukereg, ff, j, g1, wmin1;
local L_bound, U_bound, M_point;

Nobs=rows(uhat);
ukereg=zeros(Nobs,1);
Nobs=rows(uhat);
bootidx = 1+trunc(Nobs*rndus(nobs2,1,iseed));
x_boot = x[bootidx,.];
uhatboot = uhat[bootidx];

L_bound = maxc(stdc(x))*inv(rows(uhat))^(inv(cols(x)+6));
U_bound = rows(uhat)*minc(stdc(x));
M_point = L_bound+(U_bound-L_bound)*0.61803399;
{g1,wmin1} =
search(L_bound,M_point,U_bound,.05,uhatboot,x_boot,Nobs2,2);

j = 1 ;
do while j <= Nobs ;

    bootidx = 1+trunc(Nobs*rndus(nobs2,1,iseed)) ;
    x_boot = x[bootidx,.] ;
    uhatboot = uhat[bootidx] ;
    ff = ukerreg(uhatboot,x_boot,wmin1,x[j,.],&nahm,.5) ;
    ukereg[j] = ff[1,1]  ;

    j = j + 1 ;

endo ;

retp(ukereg);
endp;


proc ktest(x,uhat,Nobs2,iseed);
 local g_hat, out;
 g_hat=boot(x,uhat,Nobs2,iseed);
 out=test(uhat,g_hat,x);
retp(out);
endp;


proc test(uhat,ghat,x) ;

local A,AA,c,cc,d,n,s2,s,top,w ;

n = rows(x) ;
A = (1/n).*(x'x) ;
AA = inv(A) ;
d = uhat^2 ;
c = ghat-(x*(AA*meanc(ghat.*x))) ;
cc = c^2 ;
s = d.*cc ;
s2 = meanc(s);

top = meanc(uhat.*ghat);

w = n.*(top^2/s2) ;

retp(w) ;
endp ;

@***********************************************************************
            END OF TEST SECTION
************************************************************************@
