            PROGRAM MCARLO1
            IMPLICIT NONE
            DOUBLE PRECISION DSEED,gama1,gama2
            DOUBLE PRECISION AJGORRO,X1GORRO,X2GORRO,YGORRO
            DOUBLE PRECISION EPSX1(100),EPSX2(100),EPSY(100)
            DOUBLE PRECISION AJ(100),EPSAJ(100),BETA2,SD2,ROB2
            DOUBLE PRECISION VARZ,TRATIO,DA,P2,B1,P1,Q1,Q2
            DOUBLE PRECISION X1(100),X2(100),Y(100),Z(100)
            DOUBLE PRECISION U(100),V1(100),V2(100)
            DOUBLE PRECISION USTAR(100),V1STAR(100),V2STAR(100)
            DOUBLE PRECISION ZSSTAR(100),PI,VARR
            real mu,poidev
            INTEGER II,JJ,REP,N,I,LF,M,S
            INTEGER kn,idum
            OPEN(9,FILE='M1K2.OUT')
            OPEN(10,FILE='M1K5.OUT')
            OPEN(11,FILE='M1K8.OUT')

            mu=2.
            idum=-1
            
          REP=1000
            PI=3.1415927
           N=100
           GAMA1=2./3.
           GAMA2=SQRT(5.)/3.
           

            DSEED=2444.
            DO 10 JJ=1,REP
                 CALL GGUBS(DSEED,N,USTAR)
                 CALL GGUBS(DSEED,N,V1STAR)
                 CALL GGUBS(DSEED,N,V2STAR)
                 CALL GGUBS(DSEED,N,ZSSTAR)


              

               
               

               DO 20 I=1,N
                  U(I)=SQRT((-2)*LOG(USTAR(I)))*COS(2*ZSSTAR(I)*PI)
                  V1(I)=SQRT((-2)*LOG(V1STAR(I)))*COS(2*V2STAR(I)*PI)
                  V2(I)=SQRT((-2)*LOG(V1STAR(I)))*SIN(2*V2STAR(I)*PI)
                 Z(I)=poidev(mu,idum)
                 X1(I)=Z(I)+V1(I)
                 X2(I)=Z(I)+gama1*v1(i)+gama2*V2(I)
                     Y(I)=X1(I)+Z(I)+U(I)
20                 CONTINUE
                   DO 50 II=1,3
                   IF (II.EQ.1) kn=2
                   IF (II.EQ.2) kn=5
                   IF (II.EQ.3) kn=8



C====== Semiparametric Estimate under the alternative



               DO 200 LF=1,N
            call aknn(x2,z,kn,n,lf,x2gorro)
                    EPSX2(LF)=X2(LF)-X2GORRO

            call aknn(y,z,kn,n,lf,ygorro)
                    EPSY(LF)=Y(LF)-YGORRO
200            continue


                    CALL OLSU(EPSY,EPSX2,n,BETA2,SD2,ROB2,varz)



C======== Fitted values under the alternative

               DO 300 M=1,N
                     AJ(M)=BETA2*X2(M)
300            CONTINUE
C======== Semiparametric estimate under the composite

               DO 400 S=1,N

            call aknn(x1,z,kn,n,s,x1gorro)
                     EPSX1(S)=X1(S)-X1GORRO

            call aknn(aj,z,kn,n,s,ajgorro)
                     EPSAJ(S)=AJ(S)-AJGORRO
400             CONTINUE


            call OLSB(EPSY,EPSX1,EPSAJ,n,b1,da,p1,p2,q1,q2,VARR)
                  TRATIO=da/P2


             IF (II.EQ.1) THEN
               write(9,*) tratio
                     endif
                   IF (II.EQ.2) THEN
               write(10,*) tratio
                     endif
                   IF (II.EQ.3) THEN
               write(11,*) tratio
                     endif



50           CONTINUE
10          CONTINUE


            STOP
            END






C======== Ordinary Least Squares - Univariate and Bivariate Estimation

              SUBROUTINE OLSU(Y,X1,NOBS,BETA,SD,ROB,var)
              IMPLICIT NONE
              DOUBLE PRECISION Y(1),X1(1),BETA,SD,ROB,RS
              DOUBLE PRECISION E,EEE,RES,SUM,VAR
              DOUBLE PRECISION RSF,EF
              INTEGER J,JJ,nobs
              E=0.
              EEE=0.
              DO 238 J=1,NOBS
                  E=E+X1(J)*X1(J)
                  EEE=EEE+X1(J)*Y(J)
238           CONTINUE
              BETA=0.
              IF (E.GT.0) BETA=EEE/E
              SUM=0.
              RS=0.
              DO 239 JJ=1,NOBS
              RES=Y(JJ)-BETA*X1(JJ)
              SUM=SUM+RES*RES
              RS=RS+RES*RES*X1(jj)*X1(jj)
239           CONTINUE
              VAR=SUM/FLOAT(NOBS)
              EF=E/FLOAT(NOBS)
              RSF=RS/FLOAT(NOBS)
              SD=0.
              ROB=0.
              IF (E.GT.0) SD=SQRT(VAR/(FLOAT(NOBS)*EF))
              IF (E.GT.0) ROB=SQRT(RSF/(FLOAT(NOBS)*EF*EF))
              RETURN
              END

          SUBROUTINE OLSB(Y,X,Z,N,B1,B2,S1,S2,R1,R2,VAR)
                IMPLICIT NONE
                DOUBLE PRECISION Y(1),X(1),Z(1),B1,B2
                DOUBLE PRECISION S1,S2,R1,R2,VAR,SXX,SZZ
                DOUBLE PRECISION SXY,SZY,SRES,SR1,SR2,SR3
                DOUBLE PRECISION E1,E2,E3,ES1,ES2,ES3,FN
                DOUBLE PRECISION SXZ,DET,RES,RR1,RR2
                INTEGER N,J,K

                SXX=0.
                SXZ=0.
                SZZ=0.
                SXY=0.
                SZY=0.
                DO 100 J=1,N
                   SXX=SXX+X(J)*X(J)
                   SXZ=SXZ+X(J)*Z(J)
                   SZZ=SZZ+Z(J)*Z(J)
                   SXY=SXY+X(J)*Y(J)
                   SZY=SZY+Z(J)*Y(J)
100              CONTINUE
                 B1=0.
                 B2=0.
                 DET=SXX*SZZ-SXZ*SXZ
                 IF (DET.GT.0.) B1=(SZZ*SXY-SXZ*SZY)/DET
                 IF (DET.GT.0.) B2=(SXX*SZY-SXZ*SXY)/DET
                 SRES=0.
                 SR1=0.
                 SR2=0.
                 SR3=0.
                 DO 200 K=1,N
                     RES=Y(K)-B1*X(K)-B2*Z(K)
                     SRES=SRES+RES*RES
                     SR1=SR1+RES*RES*X(K)*X(K)
                     SR2=SR2+RES*RES*X(K)*Z(K)
                     SR3=SR3+RES*RES*Z(K)*Z(K)
200              CONTINUE
                 FN=FLOAT(N)
                 VAR=SRES/(FN-2.)
                 E1=0.
                 E2=0.
                 E3=0.
                 S1=0.
                 S2=0.
                 R1=0.
                 R2=0.
                 IF (DET.GT.0.) E1=FN*SZZ/DET
                 IF (DET.GT.0.) E2=-FN*SXZ/DET
                 IF (DET.GT.0.) E3=FN*SXX/DET
                 ES1=SR1/(FN-2.)
                 ES2=SR2/(FN-2.)
                 ES3=SR3/(FN-2.)
                 IF (DET.GT.0.) S1=SQRT(VAR*E1/FN)
                 IF (DET.GT.0.) S2=SQRT(VAR*E3/FN)
             IF (DET.GT.0.) RR1=E1*E1*ES1+E2*E2*ES3+2.*E1*E2*ES2
             IF (DET.GT.0.) RR2=E3*E3*ES3+E2*E2*ES1+2.*E2*E3*ES2
                  IF (DET.GT.0.) R1=SQRT(RR1/FN)
                  IF (DET.GT.0.) R2=SQRT(RR2/FN)
                  RETURN
                  END

C   IMSL ROUTINE NAME   - GGUBS
C
C-----------------------------------------------------------------------
C
C   COMPUTER            - IBM/SINGLE
C
C   LATEST REVISION     - JUNE 1, 1980
C
C   PURPOSE             - BASIC UNIFORM (0,1) PSEUDO-RANDOM NUMBER
C                           GENERATOR
C
C   USAGE               - CALL GGUBS (DSEED,NR,R)
C
C   ARGUMENTS    DSEED  - INPUT/OUTPUT DOUBLE PRECISION VARIABLE
C                           ASSIGNED AN INTEGER VALUE IN THE
C                           EXCLUSIVE RANGE (1.D0, 2147483647.D0).
C                           DSEED IS REPLACED BY A NEW VALUE TO BE
C                           USED IN A SUBSEQUENT CALL.
C                NR     - INPUT NUMBER OF DEVIATES TO BE GENERATED.
C                R      - OUTPUT VECTOR OF LENGTH NR CONTAINING THE
C                           PSEUDO-RANDOM UNIFORM (0,1) DEVIATES
C
C   PRECISION/HARDWARE  - SINGLE/ALL
C
C   REQD. IMSL ROUTINES - NONE REQUIRED
C
C   NOTATION            - INFORMATION ON SPECIAL NOTATION AND
C                           CONVENTIONS IS AVAILABLE IN THE MANUAL
C                           INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C   COPYRIGHT           - 1980 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C   WARRANTY            - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN
C                           APPLIED TO THIS CODE. NO OTHER WARRANTY,
C                           EXPRESSED OR IMPLIED, IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
      SUBROUTINE GGUBS (DSEED,NR,R)
C                                  SPECIFICATIONS FOR ARGUMENTS
      INTEGER            NR
      DOUBLE PRECISION   R(NR)
      DOUBLE PRECISION   DSEED
C                                  SPECIFICATIONS FOR LOCAL VARIABLES
      INTEGER            I
      DOUBLE PRECISION   D2P31M,D2P31
C                                  D2P31M=(2**31) - 1
C                                  D2P31 =(2**31)(OR AN ADJUSTED VALUE)
      DATA               D2P31M/2147483647.D0/
      DATA               D2P31/2147483648.D0/
C                                  FIRST EXECUTABLE STATEMENT
      DO 5 I=1,NR
         DSEED = DMOD(16807.D0*DSEED,D2P31M)
    5 R(I) = DSEED / D2P31
      RETURN
      END

C======== Nonparametric Kernel Estimation


               SUBROUTINE aknn(X,T1,kn,NOBS,I,XGORRO)
               implicit none
               DOUBLE PRECISIon X(1),R,C,D
               double precision t1(1)
               DOUBLE PRECISION W1,XGORRO
               integer j,nnc,nnd,kn,nobs,i
               double precision fnnc,fnnd,fkn,rro
               double precision dis(100),ro(100)
               R=T1(I)
               C=0.
               D=0.
               XGORRO=0.
                   do 301 j=1,nobs
                     dis(j)=abs(t1(j)-r)
                     ro(j)=dis(j)
301                continue
              call sort(nobs,ro)
                 rro=ro(kn)
                 nnc=0
                 nnd=0
                   do 302 j=1,nobs
                     if (dis(j).lt.rro) nnd=nnd+1
                     if (dis(j).eq.rro) nnc=nnc+1
302                continue
                  fnnc=float(nnc)
                  fnnd=float(nnd)
                  fkn=float(kn)
                    do 303 j=1,nobs
                     w1=0.
                     if (dis(j).lt.rro) w1=1./fkn
                   if (dis(j).eq.rro) then
                     if(rro.eq.0) then
                       w1=(fkn-fnnd)/(fkn*(fnnc-1.))
                     else
                       w1=(fkn-fnnd+1)/(fkn*fnnc)
                 endif
                endif
                   if(j.eq.i) w1=0.
               D=D+W1*X(j)
303            CONTINUE
                XGORRO=D
              RETURN
               END




            SUBROUTINE SORT(N,RA)
            implicit none
            double precision ra(100),rra
            integer l,i,n,j,ir
           L=N/2+1
           IR=N
10         CONTINUE
               IF(L.GT.1) THEN
                 L=L-1
                 RRA=RA(L)
              ELSE
                 RRA=RA(IR)
                RA(IR)=RA(1)
                  IR=IR-1
                  IF (IR.EQ.1) THEN
                    RA(1)=RRA
                    RETURN
                  ENDIF
              ENDIF
              I=L
              J=L+L
20              IF (J.LE.IR) THEN
                  IF (J.LT.IR) THEN
                  IF (RA(J).LT.RA(J+1)) J=J+1
                  ENDIF
                  IF (RRA.LT.RA(J)) THEN
                  RA(I)=RA(J)
                 I=J
                J=J+J
                ELSE
                J=IR+1
               ENDIF
            GO TO 20
            ENDIF
            RA(I)=RRA
          GO TO 10
          END




      FUNCTION poidev(xm,idum)
      INTEGER idum
      REAL poidev,xm,PI
      PARAMETER (PI=3.141592654)
CU    USES gammln,ran1
      REAL alxm,em,g,oldm,sq,t,y,gammln,ran1
      SAVE alxm,g,oldm,sq
      DATA oldm /-1./
      if (xm.lt.12.)then
        if (xm.ne.oldm) then
          oldm=xm
          g=exp(-xm)
        endif
        em=-1
        t=1.
2       em=em+1.
        t=t*ran1(idum)
        if (t.gt.g) goto 2
      else
        if (xm.ne.oldm) then
          oldm=xm
          sq=sqrt(2.*xm)
          alxm=log(xm)
          g=xm*alxm-gammln(xm+1.)
        endif
1       y=tan(PI*ran1(idum))
        em=sq*y+xm
        if (em.lt.0.) goto 1
        em=int(em)
        t=0.9*(1.+y**2)*exp(em*alxm-gammln(em+1.)-g)
        if (ran1(idum).gt.t) goto 1
      endif
      poidev=em
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *5=-1i.



      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *5=-1i.



      FUNCTION gammln(xx)
      REAL gammln,xx
      INTEGER j
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6)
      SAVE cof,stp
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0,
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2,
     *-.5395239384953d-5,2.5066282746310005d0/
      x=xx
      y=x
      tmp=x+5.5d0
      tmp=(x+0.5d0)*log(tmp)-tmp
      ser=1.000000000190015d0
      do 11 j=1,6
        y=y+1.d0
        ser=ser+cof(j)/y
11    continue
      gammln=tmp+log(stp*ser/x)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *5=-1i.
