c       revised nov. 20 by petra

c	This subroutine performs local linear regression of y(nx1) on X(nx1)
c	given dimension n, and bandwith h.  It calculates mhat at X2 points
c 	based on the algorithm in Fan '93--equations 2.2 - 2.4.
c	    -X2 and n2 are a vector of additional values and the length.
c	    -mhat is a vector of fitted values at the additional values. 
c       This program only calculates fitted values at the points X2.
c

	SUBROUTINE locreg4(X,y,n,X2,n2,h,mhat,denom2)
	
c	Variable Declarations

	INTEGER n, n2
	DOUBLE PRECISION X(n), y(n), X2(n2), h
	DOUBLE PRECISION  mhat(n2)

	INTEGER j, o
	DOUBLE PRECISION dist(n), w(n), K(n)
	DOUBLE PRECISION num, denom, sn1, sn2, arg, arg2, arg3
	DOUBLE PRECISION denom2(n2)

c----------------------------------------------------------------

c	Loop over each observation to calculate fitted values
c	for each local linear regression on additional values.

	DO 2 o=1,n2

c	Initialize variables to 0 for each regression.
 
	sn1 = 0.0
	sn2 = 0.0
	mhat(o) = 0.0
	DO 110 j=1,n
	  K(j) = 0.0
	  w(j) = 0.0
	  dist(j) = 0.0
	  arg = 0.0
	  dist(j) = (X2(o)-X(j))
	  arg = dist(j) / h 
	  if(abs(arg) .LE. 1.0) then
             K(j) = (15.0/16.0)*(arg**2.0-1.0)**2.0
          end if
          sn1 = sn1 + K(j) * dist(j)
	  sn2 = sn2 + K(j) * dist(j) * dist(j)

110	CONTINUE

c	------------------------------------------------------------
c	Eq 2.4: Calculate s_n,1 and s_n,2
c	Eq. 2.3: Calculate weights, w(j), for local regression at X2_o.
c	Eq. 2.2: Calculate mhat(X_o)

	num = 0.0
	denom = 1.0 / (n * n * n)
c	denom = 0.0

	DO 180 j=1,n
	  w(j) = K(j) * (sn2 - (dist(j) * sn1))
	  num = num + (w(j) * y(j))
	  denom = denom + w(j)
	  denom2(o) = denom2(o) + w(j)
180	CONTINUE
	mhat(o) = num / denom

c	End of loop for mhat(X_o)

2  	CONTINUE

	RETURN
	END