@ *************************************************************

  PROGRAM TO CARRY OUT NONPARAMETRIC TEST OF EXOGENEITY
  IN ONE-DIMENSIONAL CASE USING REAL DATA.  

  THIS PROGRAM USES A FOURIER REPRESENTATION OF THE T OPERATOR
  and AN ANALYTIC APPROXIMATION for THE ASYMPTOTIC DISTRIBUTION.

  THIS PROGRAM USES AN ORDINARY KERNEL, NOT A BOUNDARY KERNEL.
  


  MODEL:  Y = G(X) + U; E(U|W) = 0; HYPOTH:  G(X) = E(Y|X).

  24 JUNE 2004

**************************************************************@

new;
library pgraph;
clear grid, jvec, hx;
graphset;
format 8,4;

rndseed 9858289; 
tstart = date;


alph = 0.05;  @ ***** Nominal level of test ***** @
h10 = 0.10;   @ ***** Kernel bandwidth for f_xw  ***** @
h20 = 0.10;   @ ***** Kernel bandwidth for G(X) ***** @
nf = 100;  @ ***** No. of Fourier Coefficients to Use ***** @
neig = 25; @ ***** No. of Eigenvalues to Use in Computing Critical Value ***** @
csim = 100000;  @ ***** Sample size for computing critical value ***** @
jgrph = 1;  @ ***** 1 to plot results, 0 otherwise ***** @
dgnst = 0;  @ ***** 1 for diagnostic changes, 0 otherwise ***** @

sq2 = sqrt(2);
jf = seqa(1,1,nf);
n95 = 0.95*csim;
n94 = 0.94*csim;
n93 = 0.93*csim;
n90 = 0.90*csim;

output file = c:\gauss50\npiv\Etest\output.out reset;

tstart = date;


@ ***** READ DATA HERE ***** @

dataset = PATH TO DATA HERE;
 open fy=^dataset;
 n = rowsf(fy);
 yz=readr(fy,n);
 close(fy);

"Original size:  ";; n;

@ ***** SELECT VARIABLES Y, X, W HERE ***** @

 
 

@ ***** TRANSFORM TO UNIT INTERVAL ***** @
 
x = (x - minc(x))./(maxc(x) - minc(x));
w = (w - minc(w))./(maxc(w) - minc(w));


"";
"This version uses analytic approximation to asymptotic distribution";
"This version uses an ordinary kernel, not a boundary kernel."
"";

"Sample size:  ";;  n;
"No. of Eigenvalues:  ";; neig;


output off;


hx1 = h10;
hx2 = h20;


do until hx2 > 0.2;

@ ***** EVALUATE KERNEL FUNCTION for COMPUTATION UHAT and GHAT. ***** @


 dx = (x - x')/hx2;
 fxmat = kfunc(dx);
 tmp = zeros(n,1);
 fxmat = diagrv(fxmat, tmp); 
 denx = meanc(fxmat)/hx2;

 

@ ***** COMPUTE FOURIER COEFFICIENTS OF KERNEL FUNCTION and f(x,w) ***** @
@ ***** R(I,J) IS J'TH COEFFICIENT OF (1/H)*K[(X - X(I))/H] ***** @

 r = fcoef(x);
 s = fcoef(w);

 fmat = r's/n;  @ ***** Coeffs of f(x,w) ***** @

 thmat = fmat*fmat';  @ ***** Coeffs of t(x,z) ***** @;
 
 fw = sq2*sin(pi*jf*w');  @ ***** Basis functions for W ***** @
 fx = sq2*sin(pi*jf*x');  @ ***** Basis functions for x ***** @

 bphj = fmat*fw;
 bphi = (fmat*fw - thmat*fx./(denx' + 1.e-8));

ghat = sumc(y.*fxmat)./(sumc(fxmat) + 1.e-8);


@ ***** KERNEL ESTIMATION OF G(X) -- CONDITIONAL MEAN FUNCTION ***** @


 gvec = y - ghat; 

@ ***** CARRY OUT TEST BASED ON PARAMETRIC MODEL ***** @

  ww = ones(n,1)~w~(w^2);
  xx = ones(n,1)~x~(x^2);
  inv1 = inv(ww'xx)*ww';
  inv2 = invpd(xx'xx)*xx';
  invt = inv(ww'xx)*(ww'ww)*inv(ww'xx) - invpd(xx'xx);
  ivbhat = inv1*y;
  bhat = inv2*y;

  uhat = y - xx*bhat;
  su2 = (stdc(uhat))^2;
  invt = su2*invt;
  prstat = (bhat[2:3] - ivbhat[2:3])'inv(invt[2:3,2:3])*(bhat[2:3] - ivbhat[2:3]);
 
 
@ ***** COMPUTE FOURIER VERSION OF TEST STATISTIC ***** @


 tstat = gvec'(bphj'bphj)*gvec/n;
 
 uhat = gvec;
 nu = rows(uhat);
 itmp = eye(nu);
 uhmat = diagrv(itmp,(uhat^2));

 su2 = (stdc(uhat))^2;
 
 mam = su2*bphi*bphi'/n;  

 cval = 0;

 tvals = eigh(mam);
 tvals = rev(tvals);
 tvals = tvals[1:neig];

 cmat = rndn(100000,neig); 

 cmat = cmat.*cmat;
 tn = cmat*tvals;
 tn = sortc(tn,1);
 cval95 = tn[n95];
 cval94 = tn[n94];
 cval93 = tn[n93];
 cval90 = tn[n90];


 
  output on;
 "Bandwidths, Test Stat, Cvals, Hstat:  ";; hx1;; hx2;; tstat;; cval95;; cval94;; cval93;; cval90;; prstat;
 trun = ethsec(tstart,date)/6000;
 "Running time in minutes:  ";; trun;
 output off;

  hx2 = hx2 + 0.001;
  endo;

  
  

/************ END OF MAIN PROGRAM ********************************/



 proc(1) = kfunc(x);

/****************************************************************

  QUARTIC KERNEL

*****************************************************************/

 local temp;

 temp = (15/16)*((1 - x^2)^2).*(abs(x) .<= 1);

 retp(temp);
 endp;

/*****************************************************************/

 proc fcoef(x);

@ ***** COMPUTE FOURIER COEFFIENTS OF KERNEL FUNCTION. R(I,J) IS
        J'TH COEFFICIENT OF (1/H)*K[(X - X(I))/H].  X IS COLUMN
        VECTOR ***** @

local trm1, trm2, trm3, denom, abnd, bbnd,atrm, btrm;

bbnd = (x + hx1 .> 1) + (x + hx1).*(x + hx1 .<= 1);
abnd = (x - hx1).*(x - hx1 .> 0); 

denom = pi*hx1*jf';
atrm = pi*abnd*jf';
btrm = pi*bbnd*jf';
trm1 = 0.75*(1 - (x/hx1)^2).*(cos(atrm) - cos(btrm));
trm1 = trm1./denom;
trm2 = 1.5*x./(hx1*denom^2);
trm2 = trm2.*(sin(btrm) - btrm.*cos(btrm) - sin(atrm) + atrm.*cos(atrm));
trm3 = 0.75./(denom^3);
trm3 = trm3.*(2*btrm.*sin(btrm) - (btrm^2 - 2).*cos(btrm) - 2*atrm.*sin(atrm)
         + (atrm^2 - 2).*cos(atrm));
retp(sq2*(trm1 + trm2 - trm3));
endp;