        program cvmbed
        implicit double precision(a-h,o-z)
        parameter(jdmax=10,nxmax=1000)
        dimension xt(nxmax),xdat(nxmax),xlag(nxmax,jdmax),hvals(20)
        character*30 fname,outname
        common /xdata/ xt,xlag,cgcv,jd,nx,m,ipr,jw
	external rms

c----------- Read in parameters
        open(11,file='cvmbed.par',status='old')
        read(11,10) fname
        read(11,10) outname
10      format(a30)
        read(11,*) nx,jt1,ltr,ifreq,jd1,jd2,l1,l2,
     *       jtp1,jtp2,cgcv,ipr,jw
        close(11)
        csav=cgcv
        if (jw.lt.0) jw=0

c------------ read in data series; set nx=total series length
       call getx(fname,xdat,nxmax,nx,ifreq,ltr)

c---------- Set up the output files
        open(12,file=outname,status='unknown')
        write(12,110) fname,ltr
110     format(1x,'@ ',a30,'Output from CVEMBED.FOR',3x,'ltr=',i2)
        write(12,*) 'Data or transformed data scaled to Std Dev ==1'
        write(12,120) nx,jt1
120     format(1x,'max n= ',i6,' data pts, window length= ',i5)
        write(12,125) ifreq,cgcv,jw
125     format(1x, 'ifreq= ',i4, '  Cgcv= ',f5.2,' W= ',i4)
        write(12,*) 'output tp, l, d, CV_opt, h_opt, RMS_opt, r^2_opt @'
        close(12)

c------------------- Scale so SD of values to predict == 1
       sd=sdev(xdat,jt1+1,nx)
       call dscal(nx,1./sd,xdat,1)
       xbar=0.
       do 140 j=1,nx
         xbar=xbar+xdat(j)
140    continue
       xbar=xbar/dfloat(nx)
       do 150 j=1,nx
         xdat(j)=xdat(j)-xbar
150    continue

c------------------ START MAIN LOOP on prediction time, time-delay, #lags
       do 520 jtp=jtp1,jtp2
       do 510 ldelay=l1,l2
          do 500 jd=jd1,jd2

          maxback=jtp+(jd-1)*ldelay
          if(jt1.lt.maxback) goto 500

          call recontp(xt,xdat,xlag,nxmax,jdmax,nx,ldelay,jd,jt1,
     *        m,jtp)

c
c try to bracket a minimum of CV(h)
c
	  cvmin=100.
	  imax=15
	  do 400 i=1,imax
		hi=10.**( -4.0+ 0.5*dfloat(i-1))
		hvals(i)=hi
		cvi=rms(hi)
		if (cvi.lt.cvmin) then
			cvmin=cvi
			imin=i
		endif
400	  continue

c
c if bracketing failed, use endpoint as optimal bandwidth
c
	  if (imin.eq.1) then
		hcv=hvals(1)
		cvfn=cvmin	
	  elseif (imin.eq.imax) then
		hcv=hvals(imax)
		cvfn=cvmin
	  else
c
c if bracketing succeeded, find the minimum
c
	  	tol=.00001
		ax=hvals(imin-1)
		bx=hvals(imin)
		cx=hvals(imin+1)
	   	cvfn=ffmin(ax,bx,cx,rms,tol,hcv)
	  endif
c
c save CV(h_opt)
c
          cgcv=1.
          rmscv=rms(hcv)
	  r2=1.d0-rmscv**2
          cgcv=csav
	  open(12,file=outname,status='old',access='append') 
          write(12,210) jtp,ldelay,jd,cvfn,hcv,rmscv,r2
	  close(12)
210       format(1x,3i3,4(2x,f9.4))

500    continue
510    continue
520    continue
c---------------------- END MAIN LOOP
       close(12)
       stop
       end

c------------------------------------------------------------------------
c        Leave-one-out kernel regression estimates
c------------------------------------------------------------------------
      subroutine kregfhat(h,n,m,xk,y,npts,nvar,f,params,jw)

c   Computes leave-one-out kernel regression estimates at each data point,
c   for use in cross-validation.
c
c   Kernel used is a rational approximation to multivariate normal, based on
c   exp(-(x^2))=1/(1 + 0.5*x^2+ 0.3*x^6); here x is the Euclidean distance
c   between data points
c
c   Inputs:
c       h               bandwidth
c       xk(n,m)         independent variables data (column=variable)
c       y(n)            dependent variable data
c       nvar            number of indep variables
c       npts            number of data points
c   Outputs:
c       f(n)            predicted values by leave-one-out


      implicit double precision(a-h,o-z)
      parameter(zero=0.0,one=1.0)
      dimension xk(n,m),y(n),f(n)

      rh2= one/(h*h)
c----------- Loop on j is to predict at jth data point
      do 20 J=1,npts
        params=0.0
        tempd= 0.0
        tempn= 0.0

        do 10 I= 1,npts
          if (abs(i-j).le.jw) goto 10

c-------   loop over components of the vectors 
          dist=0.0
          do 5 jj=1,nvar
              dist= dist + (xk(j,jj)-xk(i,jj))**2
5         continue

c------     accumulate contribution from each data point
          x2=dist*rh2
          tempw= one/(one+ x2+ 0.5*x2**2 + 0.3*x2**3)
          tempn= tempn + tempw*y(i)
          tempd= tempd + tempw
10      continue

c------    normalize estimate
        f(j)= tempn/tempd

c------    accumulate contribution to effective # of parameters
        params=params+ (one/(one+tempd))

20    continue

      return    
      END    

c-------------------------------------------------------------------------
c                 Cross-validation function
c-------------------------------------------------------------------------
      function rms(h)
      implicit double precision(a-h,o-z)
      parameter(zero=0.0,one=1.0,nxmax=1000,jdmax=10)
      dimension xlag(nxmax,jdmax),xt(nxmax),fhat(nxmax)
      common /xdata/ xt,xlag,cgcv,jd,nx,m,ipr,jw
c
c  Compute the leave-out-2W+1 kernel estimate at each point
c

      call kregfhat(h,nxmax,jdmax,xlag,xt,m,jd,fhat,params,jw)
c
c   subroutine kregfhat(h,n,m,xk,y,npts,nvar,f,params,jw)
c       h               bandwidth
c       xk(n,m)         independent variables data (column=variable)
c       y(n)            dependent variable data
c       npts            number of data points
c       nvar            number of indep variables
c       jw              W, for leave-out-2W+1 estimate
c   Outputs:
c       f(n)            predicted values by leave-out-2W+1
c       params          effective number of parameters, tr(A(h)) 

      sse=0.
      do 2 j=1,m
        sse=sse+(fhat(j)-xt(j))**2
2     continue
      rms=sqrt(sse/dfloat(m-1))
      if (ipr.gt.0) write(*,*) rms,h,params
c
c       Compute modified cross-validation function V_c
c
      pmax=dfloat(m)/cgcv
      if (params.ge.pmax) then
         rms=10.
      else
         xm=dfloat(m)
         rms=rms*(xm-params)/(xm-cgcv*params)
      endif
      return
      end
