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

  TEST OF PARAMETRIC MODEL VS NONPARAMETRIC ALTERNATIVE WITH 
  IV ESTIMATION IN ONE-DIMENSIONAL CASE.

  THIS PROGRAM USES A FOURIER REPRESENTATION OF THE T OPERATOR
  and AN ANALYTIC APPROXIMATION for THE ASYMPTOTIC DISTRIBUTION.
  
 
  MODEL:  Y = G(X) + U; E(U|W) = 0; HYPOTH:  G(X) = B0 +X*B1

  22 JUNE 2004

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

new;
library pgraph;
clear grid, jvec, hx, oval;
graphset;
format 9,5;

rndseed 9858389; 

tstart = date;

alph = 0.05;  @ ***** Nominal level of test ***** @
hx = 0.06;   @ ***** Kernel bandwidth for test statistic ***** @
dh = 0.01;   @ ***** Increment in bandwidth ***** @
hmax = 0.06;   @ ***** Maximum bandwidth ***** @
nf = 200;  @ ***** No. of Fourier Coefficients to Use ***** @
neig = 25; @ ***** No. of Eigenvalues to Use to Compute Critical Value ***** @
iquad = 0;  @ ***** 1 for quadratic null hypothesis.  0 for linear null ***** @



 output file = c:\gauss50\npiv\yourfile.out on;  

@ ***** READ THE DATA HERE and STORE THEM AS Y, X, W.**** @



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

xx = ones(n,1)~x;
ww = ones(n,1)~w;


if iquad == 1;
 xx = ones(n,1)~x~(x^2);
 ww = ones(n,1)~w~(w^2);
 endif;

"Sample size:  ";;  n;
"No. of Fourier Coefficients: ";; nf;
"No. of Eigenvalues in Crit. Val. Computation:  ";; neig;

output off;


sq2 = sqrt(2);
jvec = seqa(1,1,nf);

z = seqa(.05,.05, 19);
rz = rows(z);
jf = seqa(1,1,nf);

tstart = date;



@ ***** 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) ***** @
 tmat = fmat*fmat';  @ ***** Fourier Coeffs of T operator ***** @

 tvals = eigh(tmat);  @ ***** Eigenvalues of T operator ***** @
 tvals = rev(tvals);
 sel = (tvals .> 0);
 

 fw = sq2*sin(pi*jf*w');  @ ***** Basis functions for W ***** @
 bphi = fmat*fw;

@ ***** EVALUATE KERNEL FUNCTION for COMPUTATION OF UHAT ***** @

 dw = (w - w')/hx;
 amat = kfunc(dw);


@ ***** IV Estimation of Parameter ***** @

  bhat = inv(ww'xx)*ww'y;
  gvec = y - xx*bhat; 


@ ***** COMPUTE FOURIER VERSION OF TEST STATISTIC ***** @


 tstat = gvec'(bphi'bphi)*gvec/n;
 
@ ***** Generate Pseudo U's for Estimation of Asymp. Distr. Func. ***** @

 qrfunc = (amat*gvec)./(sumc(amat) + 1e-8);
 uhat = gvec - qrfunc;

 invwx = inv(ww'xx);
 mmat = eye(n) - xx*invwx*ww';
 
 su2 = (stdc(uhat))^2;
 
 mam = bphi*mmat*mmat'bphi';
 tvals = eigh(mam);
 

 tvals = rev(tvals);
 tvals = tvals[1:neig];
 cmat = rndn(1000,neig); 
 cmat = cmat.*cmat;
 tn = su2*cmat*tvals/n;

 tn = sortc(tn,1);
 nalph = ceil(1000*(1 - alph));
 cval = tn[nalph];
 
 output on;
 "Bandwidth and Statistic, Crit. Val.:  ";; hx;; tstat;; cval;
 trun = ethsec(tstart,date)/6000;
 "Running time in minutes:  ";; trun;
 output off;



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


 proc(1) = kfunc(x);

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

  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 + hx .> 1) + (x + hx).*(x + hx .<= 1);
abnd = (x - hx).*(x - hx .> 0); 
@x~abnd~bbnd; @
denom = pi*hx*jf';
atrm = pi*abnd*jf';
btrm = pi*bbnd*jf';
trm1 = 0.75*(1 - (x/hx)^2).*(cos(atrm) - cos(btrm));
trm1 = trm1./denom;
trm2 = 1.5*x./(hx*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;