kern_reg.f90

The subroutines glkern.f and lokern.f use an efficient and fast algorithm for automatically adapti

  ii=ii+1
  IF(iil == 0) iil=i
  171 IF(tt(j) < t(i)) THEN
    j=j+1
    IF(j <= m) GO TO 171
  END IF
  wn(ii,3)=x(i)-y(j) + (y(j)-y(j-1))*(tt(j)-t(i))/(tt(j)-tt(j-1))
END DO
IF(iil == 0 .OR. ii-iil < 10) THEN
  CALL resest(t(il:),wn(1:,3),nn,wn(1:,4),snr,rvar)
ELSE
  CALL resest(t(iil:), wn(1:,3), ii, wn(1:,4), snr, rvar)
END IF
q=sig/rvar
IF(q <= 2.) RETURN
IF(q > 5. .AND. r2 > .95) rvar=rvar*.5
sig=rvar
CALL coff(wn(1:,4), n, sig)
GO TO 100
END SUBROUTINE glkern



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

! N.B. Arguments WN, W1 & WM have been removed.
 
!-------------------------------------------------------------------*
!--------------------------------------------------------------------
!    SHORT-VERSION: JANUARY 1997

!    PURPOSE:

!    GENERAL SUBROUTINE FOR KERNEL SMOOTHING:
!    COMPUTATION OF ITERATIVE PLUG-IN ALGORITHM FOR LOCAL 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 LOCAL PLUG-IN BANDWIDTH ARRAY IS GIVEN BY BAN(1),...,BAN(M)
!-------------------------------------------------------------------*

!  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
!  INPUT    M            LENGTH OF TT

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

!  INPUT    NUE          ORDER OF DERIVATIVE (0-4)
!                        ****** 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 LOCAL BANDWIDTH
!                        <>0 USING LOCAL INPUT BANDWIDTH-ARRAY IN BAN
!                        ****** 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 BAN(M)       LOCAL PLUG-IN BANDWIDTH ARRAY
!                        ****** WILL BE SET IN SUBROUTINE FOR ISMO=0

! WORK     WN(0:N,5)     WORK ARRAY FOR KERNEL SMOOTHING ROUTINE
!                        ****** WILL BE SET IN SUBROUTINE
! WORK     W1(M1,3)      WORK ARRAY FOR INTEGRAL APPROXIMATION
!                        ****** WILL BE SET IN SUBROUTINE
! WORK     WM(M)         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)        :: 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)  :: ban(:)
REAL (dp), INTENT(OUT)     :: y(:)

! Local variables

REAL (dp)  :: wn(0:n,5), w1(m1,3), wm(m)
INTEGER    :: i, ii, iil, il, inputs, iprint, isort, it, itende, itt, iu,  &
              j, kk, kk2, kordv, nn, nuev, nyg, nyl
REAL (dp)  :: alpha, b, b2, bmax, bmin, bres, bs, bvar, const, dist,  &
              ex, exs, exsvi, fac, g1, g2, q, r2, rvar, s0, sn, snr, ssi, &
              tll, tuu, vi, wstep, xh, xi, xmy2, xxh
!-
!-------- 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 > 0
!-------- 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(*, *) 'lokern: Order of derivative not allowed'
  STOP
END IF
IF(nue > 2 .AND. ismo == 0) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Order of derivative not allowed'
  STOP
END IF
IF(n <= 2) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Number of data too small'
  STOP
END IF
IF(m < 1) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: No output points'
  STOP
END IF
IF(m1 < 10) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Variable M1 is choosen too small'
  STOP
END IF
kk=(kord-nue)/2
IF(2*kk+nue /= kord) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
END IF
IF(kord > 4 .AND. ismo == 0) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
END IF
IF(kord > 6 .OR. kord <= nue) THEN
  IF(iprint == 0) WRITE(*, *) 'lokern: Kernel order not allowed, set to ',nue+2
  kord=nue+2
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 230
ELSE
  IF(ismo /= 0) GO TO 230
END IF
!-
!-------- 3. COMPUTATION OF MINIMAL, MAXIMAL ALLOWED GLOBAL 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(*, *)   &
    'lokern: Extrapolation at both boundaries not optimized'
IF(tt(1) < s0 .AND. tt(m) <= sn .AND. iprint < 0) WRITE(*, *)   &
    'lokern: Extrapolation at left boundary not optimized'
IF(tt(1) >= s0 .AND. tt(m) > sn .AND. iprint < 0) WRITE(*, *)   &
    'lokern: 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(s0,tl)
tu=MIN(sn,tu)
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,.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(1:,3),n,wn(1:,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 230
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/(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 GLOBAL PLUG-IN BANDWIDTH
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) GO TO 180
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
  ii=ii+1
  IF(iil == 0) iil=i
  171 IF(tt(j) < t(i)) THEN
    j=j+1
    IF(j <= m) GO TO 171
  END IF
  wn(ii,3)=x(i)-y(j) + (y(j)-y(j-1))*(tt(j)-t(i))/(tt(j)-tt(j-1))
END DO
CALL resest(t(iil:), wn(1:,3), ii, wn(1:,4), snr, rvar)
q=sig/rvar
CALL coff(wn(1:,4), n, sig)
IF(q <= 2.) GO TO 180
IF(q > 5. .AND. r2 > .95) rvar=rvar*.5
sig=rvar
CALL coff(wn(1:,4), n, sig)
GO TO 100
!-
!-------- 18. LOCAL INITIALIZATIONS
180 bvar=b
nuev=0
kordv=2
nyl=1
!-
!-------- 19. COMPUTE INNER BANDWIDTHS
g1=0.86*(1.0 + DBLE(kord-nue-2)*.05)*b
g1=g1*DBLE(n)**(4./DBLE(2*kord+1)/(2*kord+5))
g1=MAX(g1, bmin/DBLE(kord-1)*DBLE(kord+1))
g1=MIN(g1, bmax)

g2=1.4*(1.0 + DBLE(kord-nue-2)*0.05)*b
g2=g2*DBLE(n)**(2./DBLE(2*kord+1)/DBLE(2*kord+3))
g2=MAX(g2, bmin)
g2=MIN(g2, bmax)
!-
!-------- 20. ESTIMATE/COMPUTE INTEGRAL CONSTANT VI LOCALLY
DO  i=1,n
  wn(i,4)=DBLE(n)*wn(i,4)*(s(i)-s(i-1))
END DO
DO  j=1,m
  ban(j)=bvar
  wm(j)=tt(j)
  IF(tt(j) < s(0)+g1) THEN
    dist=((tt(j)-g1-s(0))/g1)**2
    ban(j)=bvar*(1.0 + dist)
    ban(j)=MIN(ban(j), bmax)
    wm(j)=tt(j) + 0.5*dist*g1
  ELSE IF(tt(j) > s(n)-g1) THEN
    dist=((tt(j)-s(n)+g1)/g1)**2
    ban(j)=bvar*(1.0 + dist)
    ban(j)=MIN(ban(j), bmax)
    wm(j)=tt(j) - .5*dist*g1
  END IF
END DO
CALL kernel(t, wn(1:,4), n, bvar, nuev, kordv, nyl, s, wm, m, ban)
!-
!-------- 21. ESTIMATION OF KORDTH DERIVATIVE LOCALLY
wstep=(tt(m)-tt(1))/(m1-2)
DO  j=2,m1
  w1(j,2)=tt(1) + (j-2)*wstep
  w1(j,1)=tt(1) + (j-1.5)*wstep
END DO
w1(1,1)=tt(1) + 0.5*wstep