MODULE Kernel_Regression

! General Remarks
! The subroutines glkern.f and lokern.f use an efficient and fast algorithm for
! automatically adaptive nonparametric regression estimation with a kernel method.
! Roughly speaking, the method performs a local averaging of the observations when
! estimating the regression function.  Analogously, one can estimate derivatives
! of small order of the regression function.

! Crucial for the kernel regression estimation used here is the choice of global
! or local bandwidths. Too small ones will lead to a wiggly curve, too large ones
! will smooth away important details.

! The subroutine glkern.f calculates an estimator of the regression function or
! derivatives of the regression function with an automatically chosen global
! plugin bandwidth.  The subroutine lokern.f calculates such an estimator with an
! automatically chosen bandwidth function.  It is also possible to use global and
! local bandwidths, respectively, which are specified by the user.
! The main idea of the plugin method is to estimate the optimal bandwidths by
! estimating the asymptotically optimal mean (integrated) squared error optimal
! bandwidths.  Therefore, one has to estimate the variance for homoscedastic error
! variables and a functional of a smooth variance function for heteroscedastic
! error variables, respectively.  Also, one has to estimate an integral functional
! of the squared k-th derivative of the regression function (k=KORD) for the
! global bandwidth and the squared k-th derivative itself for the local
! bandwidths.  Here, a further kernel estimator for this derivative is used with a
! bandwidth which is adapted iteratively to the regression function.
! A convolution form of the kernel estimator for the regression function and its
! derivatives is used.  Thereby, one can adapt the S-array (which is
! S(I)=(T(I)+T(I+1))/2 in the standard convolution form) to be a smoothed grid
! more suitable for random design, see Herrmann (1996). Using this estimator leads
! to an asymptotically minimax efficient estimator for fixed and random design.
! Polynomial kernels and boundary kernels are used with a fast and stable updating
! algorithm for kernel regression estimation.

! More details can be found in the papers referred to in the references.

! References
! On the global iterative plugin bandwidth estimator:
! T. Gasser, A. Kneip, and W. K鰄ler (1991). A flexible and fast method for
! automatic smoothing. Journal of the American Statistical Association, 86,
! 643-652.
! On the local plugin bandwidth estimator:
! M. Brockmann, T. Gasser, and E. Herrmann (1993). Locally adaptive bandwidth
! choice for kernel regression estimators. Journal of the American Statistical
! Association, 88, 1302-1309.
! On nonparametric variance estimation:
! T. Gasser, L. Sroka, and C. Jennen-Steinmetz (1986). Residual and residual
! pattern in nonlinear regression. Biometrika, 73, 625-633.
! On adapting heteroscedasticity:
! E. Herrmann (1997). Local bandwidth choice in kernel regression estimation.
! Journal of Graphical and Computational Statistics, 6, 35-54.
! On the fast algorithm for kernel regression estimator:
! T. Gasser and A. Kneip (1989) discussion of Buja, A., Hastie, T. and Tibshirani,
! R.: Linear smoothers and additive models, The Annals of Statistics, 17, 532-535.
!
! B. Seifert, M. Brockmann, J. Engel, and T. Gasser (1994). Fast algorithms for
! nonparametric curve estimation. J. Computational and Graphical Statistics 3,
! 192-213.
! On the special kernel estimator for random design:
! E. Herrmann (1996). On the convolution type kernel regression estimator.
! Preprint 1833, FB Mathematik, Technische Universit鋞 Darmstadt (available from
! the preprint server of the mathematical department of the technical university
! of Darmstadt )

! Comments to Eva Herrmann: eherrmann@mathematik.tu-darmstadt.de

! Last update (Fortran 77): 10-December-98 / mz

! Code converted using TO_F90 by Alan Miller
! Date: 2000-10-04  Time: 16:41:52

IMPLICIT NONE

INTEGER, PARAMETER  :: dp = SELECTED_REAL_KIND(12, 60)


CONTAINS


SUBROUTINE glkern(t, x, n, tt, m, ihom, nue, kord, irnd,  &
                  ismo, m1, tl, tu, s, sig, b, y)

! N.B. Arguments WN & W1 have been removed.
 
!------------------------------------------------------------------*
!   SHORT-VERSION: OCT 1996

!   PURPOSE:

!   GENERAL SUBROUTINE FOR KERNEL SMOOTHING:
!   COMPUTATION OF ITERATIVE PLUG-IN ALGORITHM FOR GLOBAL BANDWIDTH
!   SELECTION FOR KERNELS WITH (NUE,KORD) = (0,2),(0,4),(1,3) OR (2,4).

!------------------------------------------------------------------*
!   THE RAW DATA SHOULD BE GIVEN BY THE POINTS
!   (T(1),X(1)),...,(T(N),X(N))

!   THE RESULTING ESTIMATOR OF THE NUE-TH DERIVATIVE OF THE
!   REGRESSION CURVE IS GIVEN THROUGH THE POINTS
!   (TT(1),Y(1)),...,(TT(M),Y(M))

!   THE PLUG-IN BANDWIDTH IS GIVEN BY B
!------------------------------------------------------------------*

!  PARAMETERS :

!  INPUT    T(N)         INPUT GRID (T(1)<T(2)<...<T(N))
!  INPUT    X(N)         DATA
!  INPUT    N            LENGTH OF X

!  INPUT    TT(M)        OUTPUT GRID, SHOULD BE ORDERED
!  INPUT    M            LENGTH OF TT

!  INPUT    IHOM         HOMOSKEDASTICY OF VARIANCE
!                        0: HOMOSKEDASTIC ERROR VARIABLES,
!                        <> 0: IF THE VARIANCE SHOULD ESTIMATED AS
!                              SMOOTH FUNCTION.
!                        ****** DEFAULT VALUE: IHOM=0

!  INPUT    NUE          ORDER OF DERIVATIVE (0-4) OF THE REGRESSION
!                        FUNCTION WHICH SHALL BE ESTIMATED
!                        ****** DEFAULT VALUE: NUE=0

!  INPUT    KORD         ORDER OF KERNEL (<=6), FOR ISMO=0  ONLY
!                        NUE=0, KORD=2 OR KORD=4
!                        OR NUE=1, KORD=3 OR NUE=2, KORD=4 ARE ALLOWED
!                        ****** DEFAULT VALUE: KORD=NUE+2

!  INPUT    IRND         0: IF RANDOM GRID POINTS T MAY OCCUR
!                        <>0 ELSE  (ONLY NECESSARY IF S SHOULD BE
!                            COMPUTED)
!                        ****** DEFAULT VALUE IRND=0

!  INPUT    ISMO         0:ESTIMATING THE OPTIMAL GLOBAL BANDWIDTH
!                        <>0 USING GLOBAL INPUT BANDWIDTH B
!                        ****** DEFAULT VALUE ISMO=0

!  INPUT    M1           >=10, LENGTH OF W1, LARGE VALUES WILL INCREASE
!                        THE ACCURACY OF THE INTEGRAL APPROXIMATION
!                        ****** DEFAULT VALUE: M1=400

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! IN/OUTPUT TL/TU        LOWER/UPPER BOUND FOR INTEGRAL APPROXIMATION
!                        AND VARIANCE ESTIMATION (IF SIG=0 AND IHOM=0),
!                        IF TU<=TL, [TL,TU] ARE COMPUTED AS ABOUT
!                        THE 87% MIDDLE PART OF [T(1),T(N)]
!                        ****** DEFAULT VALUES: TL=1.0, TU=0.0

! IN/OUTPUT S(0:N)       IF S(N)<=S(0) THIS ARRAY IS COMPUTED AS
!                        MIDPOINTS OF T, FOR NON-RANDOM DESIGN AND AS
!                        SMOOTHED QUANTILES FOR RANDOM DESIGN
!                        ****** DEFAULT VALUES: S(0)=1.0, S(N)=0.0
!                               AND THE OTHER S(I) CAN BE UNDEFINED

! IN/OUTPUT SIG          RESIDUAL VARIANCE, ESTIMATED FOR SIG=0 OR
!                        IHOM<>0, ELSE GIVEN BY INPUT
!                        ****** DEFAULT VALUE: SIG=-1.0

! IN/OUTPUT B            GLOBAL PLUG-IN BANDWIDTH
!                        ****** B CAN BE UNDIFINED IF ISMO=0


! WORK     WN(0:N,5)     WORK ARRAY FOR KERNEL SMOOTHING ROUTINE
!                        OR NUE=1, KORD=3 OR NUE=2, KORD=4 ARE ALLOWED
!                        ****** WILL BE SET IN SUBROUTINE
! WORK     W1(M1,3)      WORK ARRAY FOR INTEGRAL APPROXIMATION
!                        ****** WILL BE SET IN SUBROUTINE

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! OUTPUT   Y(M)          KERNEL ESTIMATE WITH BOP (=B0 FOR ISMO<>0)
!                        ****** WILL BE SET IN SUBROUTINE
!-----------------------------------------------------------------------
!  USED SUBROUTINES: COFF, RESEST, KERNEL WITH FURTHER SUBROUTINES
!                    WHICH ARE CONTAINED IN THE FILE subs.f
!-----------------------------------------------------------------------

REAL (dp), INTENT(IN)      :: t(:)
REAL (dp), INTENT(IN)      :: x(:)
INTEGER, INTENT(IN)        :: n
REAL (dp), INTENT(IN)      :: tt(:)
INTEGER, INTENT(IN)        :: m
INTEGER, INTENT(IN)        :: ihom
INTEGER, INTENT(IN)        :: nue
INTEGER, INTENT(IN OUT)    :: kord
INTEGER, INTENT(IN)        :: irnd
INTEGER, INTENT(IN OUT)    :: ismo
INTEGER, INTENT(IN)        :: m1
REAL (dp), INTENT(IN OUT)  :: tl
REAL (dp), INTENT(IN OUT)  :: tu
REAL (dp), INTENT(IN OUT)  :: s(0:)
REAL (dp), INTENT(IN OUT)  :: sig
REAL (dp), INTENT(IN OUT)  :: b
REAL (dp), INTENT(OUT)     :: y(:)

! Local variables

REAL (dp)  :: w1(m1,3), wn(0:n,5)
INTEGER    :: i, ii, iil, il, inputs, iprint, isort, it, itende, itt, iu,  &
              j, kk, kk2, nn, nyg
REAL (dp)  :: alpha, b2, bmax, bmin, bres, bs, const, ex, exs, exsvi, fac,  &
              q, r2, rvar, s0, sn, snr, ssi, tll, tuu, vi, xi, xmy2
!-
!-------- 1. INITIALISATIONS AND SOME ERROR-CHECKS
REAL (dp), SAVE  :: bias(2,0:2) = RESHAPE(  &
                    (/ 0.2, 0.04762, 0.4286, 0.1515, 1.33, 0.6293 /), (/ 2,3 /))
REAL (dp), SAVE  :: vark(2,0:2) = RESHAPE(  &
                    (/ 0.6, 1.250, 2.143, 11.93, 35.0, 381.6 /), (/ 2,3 /))
REAL (dp), SAVE  :: fak2(2:4) = (/ 4., 36., 576. /)

nyg=0
inputs=0
!-------- IF NO ERRORS SHOULD BE WRITTEN ON STANDARD OUTPUT, SET IPRINT=1
!-------- IF ERRORS AND VERY DETAILED WARNINGS SHOULD BE WRITTEN ON
!--------           STANDARD OUTPUT, SET IPRINT < 0
iprint=0

IF(nue > 4 .OR. nue < 0) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Order of derivative not allowed'
  STOP
END IF
IF(nue > 2 .AND. ismo == 0) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Order of derivative not allowed'
  STOP
END IF
IF(n <= 2) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Number of data too small'
  STOP
END IF
IF(m < 1) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: No output points'
  STOP
END IF
IF(m1 < 10) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Variable M1 is choosen too small'
  STOP
END IF

kk=(kord-nue)/2
IF(2*kk+nue /= kord) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
END IF
IF(kord > 4 .AND. ismo == 0) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
END IF
IF(kord > 6 .OR. kord <= nue) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
END IF
IF(ismo /= 0 .AND. b <= 0) THEN
  IF(iprint == 0) WRITE(*, *) 'glkern: Plug-in bandwidth is used'
  ismo=0
END IF
rvar=sig
!-
!-------- 2. COMPUTATION OF S-SEQUENCE
s0=1.5*t(1) - 0.5*t(2)
sn=1.5*t(n) - 0.5*t(n-1)
IF(s(n) <= s(0)) THEN
  inputs=1
  DO  i=1,n-1
    s(i)=.5*(t(i)+t(i+1))
  END DO
  s(0)=s0
  s(n)=sn
  IF(ismo /= 0 .AND. irnd /= 0) GO TO 160
ELSE
  IF(ismo /= 0) GO TO 160
END IF
!-
!-------- 3. COMPUTATION OF MINIMAL, MAXIMAL ALLOWED BANDWIDTH
bmax=(sn-s0)*.5
bmin=(sn-s0)/DBLE(n)*DBLE(kord-1)*.6
!-
!-------- 4. WARNINGS IF TT-GRID LARGER THAN T-GRID
IF(tt(1) < s0 .AND. tt(m) > sn .AND. iprint < 0) WRITE(*, *)   &
    'glkern: Extrapolation at both boundaries not optimized'
IF(tt(1) < s0 .AND. tt(m) <= sn .AND. iprint < 0) WRITE(*, *)   &
    'glkern: Extrapolation at left boundary not optimized'
IF(tt(1) >= s0 .AND. tt(m) > sn .AND. iprint < 0) WRITE(*, *)   &
    'glkern: Extrapolation at right boundary not optimized'

!-
!-------- 5. COMPUTE TL,TU AND THEIR T-GRID AS INNER PART FOR
!            INTEGRAL APPROXIMATION IN THE ITERATIONS
itt=0
51 IF (tu <= tl) THEN
  tl=.933*s0 +.067*sn
  tu=.067*s0 +.933*sn
  itt=itt+1
END IF
tl=MAX(tl,s0)
tu=MIN(tu,sn)
il=1
iu=n
wn(1,1)=0.0
wn(n,1)=0.0
DO  i=1,n
  IF(t(i) <= tl .OR. t(i) >= tu) wn(i,1)=0.0
  IF(t(i) > tl .AND. t(i) < tu) wn(i,1)=1.0
  IF(t(i) < tl) il=i+1
  IF(t(i) <= tu) iu=i
END DO
nn=iu-il+1
IF(nn == 0 .AND. itt == 0) THEN
  tu=tl-1.0
  GO TO 51
END IF
IF(nn == 0 .AND. itt == 1) THEN
  tu=sn
  tl=s0
  GO TO 51
END IF
!-
!-------- 6. COMPUTE T-GRID FOR INTEGRAL APPROXIMATION
DO  i=1,m1
  w1(i,2)=1.0
  w1(i,1)=tl+(tu-tl)*DBLE(i-1)/DBLE(m1-1)
END DO
!-
!-------- 7. CALCULATION OF WEIGHT FUNCTION
alpha=1.d0/DBLE(13)
DO  i=il,iu
  xi=(t(i) - tl)/alpha/(tu-tl)
  IF(xi > 1) GO TO 71
  wn(i,1)=(10.0 - 15*xi + 6*xi*xi)*xi*xi*xi
END DO
71 DO  i=iu,il,-1
  xi=(tu-t(i))/alpha/(tu-tl)
  IF(xi > 1) GO TO 73
  wn(i,1)=(10.0 - 15*xi + 6*xi*xi)*xi*xi*xi
END DO
73 DO  i=1,m1
  xi=(w1(i,1)-tl)/alpha/(tu-tl)
  IF(xi > 1) GO TO 75
  w1(i,2)=(10.0 - 15*xi + 6*xi*xi)*xi*xi*xi
END DO
75 DO  i=m1,1,-1
  xi=(tu-w1(i,1))/alpha/(tu-tl)
  IF(xi > 1) GO TO 77
  w1(i,2)=(10.0 - 15*xi + 6*xi*xi)*xi*xi*xi
END DO
!-
!-------- 8. COMPUTE CONSTANTS FOR ITERATION
77 ex=1./DBLE(kord+kord+1)
kk2=(kord-nue)
kk=kk2/2
!-
!-------- 9. ESTIMATING VARIANCE AND SMOOTHED PSEUDORESIDUALS
IF(sig <= .0 .AND. ihom == 0) CALL resest(t(il:),x(il:),nn,wn(il:,2),r2,sig)
IF(ihom /= 0) THEN
  CALL resest(t,x,n,wn(1:,2),snr,sig)
  bres=MAX(bmin, 0.2*nn**(-.2)*(s(iu)-s(il-1)))
  DO  i=1,n
    wn(i,3)=t(i)
    wn(i,2)=wn(i,2)*wn(i,2)
  END DO
  CALL kernel(t,wn(1:,2),n,bres,0,kk2,nyg,s,wn(il:,3),nn,wn(il:,4))
ELSE
  CALL coff(wn(1:,4),n,sig)
END IF
!-
!-------- 10. ESTIMATE/COMPUTE INTEGRAL CONSTANT
100 vi=0.
DO  i=il,iu
  vi=vi + wn(i,1)*n*(s(i)-s(i-1))**2*wn(i,4)
END DO
!-
!-------- 11. REFINEMENT OF S-SEQUENCE FOR RANDOM DESIGN
IF(inputs == 1 .AND. irnd == 0) THEN
  DO  i=0,n
    wn(i,5)=DBLE(i)/DBLE(n+1)
    wn(i,2)=(DBLE(i)+.5)/DBLE(n+1)
    wn(i,3)=wn(i,2)
  END DO
  exs=-DBLE(3*kord+1)/DBLE(6*kord+3)
  exsvi=DBLE(kord)/DBLE(6*kord+3)
  bs=0.1*(vi/(sn-s0)**2)**exsvi*n**exs
  CALL kernel(wn(1:,5),t,n,bs,0,2,nyg,wn(0:,3),wn(0:,2),n+1,s(0:))
  111 isort=0
  vi=0.0
  DO  i=1,n
    vi=vi + wn(i,1)*n*(s(i)-s(i-1))**2*wn(i,4)
    IF(s(i) < s(i-1)) THEN
      ssi=s(i-1)
      s(i-1)=s(i)
      s(i)=ssi
      isort=1
    END IF
  END DO
  IF(isort == 1) GO TO 111
  IF(ismo /= 0) GO TO 160
END IF
b=bmin*2.
!-
!-------- 12. COMPUTE INFLATION CONSTANT AND EXPONENT AND LOOP OF ITERATIONS
const=DBLE(2*nue+1)*fak2(kord)*vark(kk,nue)*vi  &
      /(DBLE(2*kord-2*nue)*bias(kk,nue)**2*DBLE(n))
fac=1.1*(1.+(nue/10.)+0.05*(kord-nue-2.)) *n**(2./DBLE((2*kord+1)*(2*kord+3)))
itende=1 + 2*kord + kord*(2*kord+1)

DO  it=1,itende
!-
!-------- 13. ESTIMATE DERIVATIVE OF ORDER KORD IN ITERATIONS
  b2=b*fac
  b2=MAX(b2,bmin/DBLE(kord-1)*DBLE(kord+1))
  b2=MIN(b2,bmax)
  CALL kernel(t,x,n,b2,kord,kord+2,nyg,s,w1(1:,1),m1,w1(1:,3))
!-
!-------- 14. ESTIMATE INTEGRALFUNCTIONAL IN ITERATIONS
  xmy2=.75*(w1(1,2)*w1(1,3)*w1(1,3) + w1(m1,2)*w1(m1,3)*w1(m1,3))
  DO  i=2,m1-1
    xmy2=xmy2 + w1(i,2)*w1(i,3)*w1(i,3)
  END DO
  xmy2=xmy2*(tu-tl)/m1
!-
!-------- 15. FINISH OF ITERATIONS
  b=(const/xmy2)**ex
  b=MAX(bmin,b)
  b=MIN(bmax,b)
END DO
!-
!-------- 16  COMPUTE SMOOTHED FUNCTION WITH PLUG-IN BANDWIDTH
160 CALL kernel(t, x, n, b, nue, kord, nyg, s, tt, m, y)
!-
!-------- 17. VARIANCE CHECK
IF(ihom /= 0) sig=rvar
IF(rvar == sig .OR. r2 < .88 .OR. ihom /= 0 .OR. nue > 0) RETURN
ii=0
iil=0
j=2
tll=MAX(tl, tt(1))
tuu=MIN(tu, tt(m))
DO  i=il,iu
  IF(t(i) < tll .OR. t(i) > tuu) CYCLE