            PROGRAM TABLE2B1
               IMPLICIT NONE
               DOUBLE PRECISION X(23972),Y(23972)
               DOUBLE PRECISION Z1(23972),Z2(23972)
               DOUBLE PRECISION Z3(23972)
               DOUBLE PRECISION Z4(23972)
               DOUBLE PRECISION XT(23972)
               DOUBLE PRECISION LX(23972),L2(23972)
               DOUBLE PRECISION EX(23972),EY(23972)
               DOUBLE PRECISION EXT(23972),MX(23972)
               DOUBLE PRECISION ELX(23972),ELX2(23972)
               DOUBLE PRECISION EAJ(23972)
               DOUBLE PRECISION AJUST(23972)
               DOUBLE PRECISION LXG,LX2G,AJGORRO
               DOUBLE PRECISION XG,XTG,YG
               DOUBLE PRECISION B1,B2,S1,S2,R1,R2,VV
               DOUBLE PRECISION G1,G2,G3,P1,P2,P3
               DOUBLE PRECISION Q1,Q2,Q3,V
               INTEGER N,J,K

               OPEN (5,FILE='y.dat')
               OPEN (6,FILE='x.dat')
               OPEN (7,FILE='z1.dat')
               OPEN (8,FILE='z2.dat')
               OPEN (9,FILE='z3.dat')
               OPEN (10,FILE='z4.dat')
               OPEN (16,FILE='table2b1.out')
               N=23972
               K=160
               DO 100 J=1,N
                   READ(5,*) Y(J)
                   READ(6,*) MX(J)
                   READ(7,*) Z1(J)
                   READ(8,*) Z2(J)
                   READ(9,*) Z3(J)
                   READ(10,*) Z4(J)

                                   
                   X(J)=LOG(MX(J))
                   XT(J)=LOG(MX(J))*LOG(MX(J))
                   LX(J)=MX(J)
                   L2(J)=1./MX(J)
 100           CONTINUE


C====== Semiparametric estimate under the alternative


               DO 200 J=1,N
       CALL A(Y,LX,L2,x,xt,Z1,z2,z3,z4,J,YG,lxg,lx2g,xg,xtg,k,n)
                    EX(J)=X(J)-XG
                    EXT(J)=XT(J)-XTG
                    EY(J)=Y(J)-YG
                     ELX(J)=LX(J)-LXG
                     ELX2(J)=L2(J)-LX2G

200            continue
                    CALL OLSB(EY,EX,EXt,N,B1,B2,S1,S2,R1,R2,VV)
                    WRITE(16,*)
                    WRITE(16,*) B1,B2
                    WRITE(16,*) S1,S2
                    WRITE(16,*) R1,R2
                    WRITE(16,*) B1/S1,B2/S2
                    WRITE(16,*) B1/R1,B2/R2

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

               DO 300 J=1,N
                     AJUST(J)=B1*X(J)+B2*XT(J)
300            CONTINUE

C======== Semiparametric estimates under the composite

               DO 400 J=1,N

       CALL Aknn(ajust,Z1,z2,z3,z4,j,ajgorro,k,n)
                     EAJ(J)=AJUST(J)-AJGORRO
400            CONTINUE

         CALL OLST(EY,ELX,ELX2,EAJ,N,G1,G2,G3,P1,P2,P3,Q1,Q2,Q3,V)
               WRITE(16,*)
               WRITE(16,*) G1,G2,G3
               WRITE(16,*) P1,P2,P3
               WRITE(16,*) Q1,Q2,Q3
               WRITE(16,*) G3/P3
               WRITE(16,*) G3/Q3
               WRITE(16,*)
               STOP
               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


        SUBROUTINE OLST(Y,X,Z,W,N,B1,B2,B3,S1,S2,S3,R1,R2,R3,VAR)
            IMPLICIT NONE
            DOUBLE PRECISION Y(1),X(1),Z(1),W(1),B1,B2,B3
            DOUBLE PRECISION S1,S2,S3,R1,R2,R3,VAR
            DOUBLE PRECISION SXX,SXZ,SXW,SZZ,SZW,SWW
            DOUBLE PRECISION SXY,SZY,SWY,DET,FN,RES
            DOUBLE PRECISION I1,I2,I3,I4,I5,I6,E1,E2
            DOUBLE PRECISION E3,E4,E5,E6,ES1,ES2,ES3,ES4
            DOUBLE PRECISION ES5,ES6,SRES,SR1,SR2,SR3,SR4
            DOUBLE PRECISION SR5,SR6,AAA,BBB,RR1,RR2,RR3
            DOUBLE PRECISION M1,M2,M3,M4,M5,M6
            INTEGER N,J,K

             SXX=0.
             SXZ=0.
             SXW=0.
             SZZ=0.
             SZW=0.
             SWW=0.
             SXY=0.
             SZY=0.
             SWY=0.
              DO 100 J=1,N
              SXX=SXX+X(J)*X(J)
              SXZ=SXZ+X(J)*Z(J)
              SXW=SXW+X(J)*W(J)
              SZZ=SZZ+Z(J)*Z(J)
              SZW=SZW+Z(J)*W(J)
              SWW=SWW+W(J)*W(J)
              SXY=SXY+X(J)*Y(J)
              SZY=SZY+Z(J)*Y(J)
              SWY=SWY+W(J)*Y(J)
100           CONTINUE
              AAA=SXX*SZZ*SWW+2.*SXZ*SZW*SXW
              BBB=SXW*SXW*SZZ+SXZ*SXZ*SWW+SZW*SZW*SXX
              DET=AAA-BBB
              I1=SZZ*SWW-SZW*SZW
              I2=SZW*SXW-SXZ*SWW
              I3=SXZ*SZW-SZZ*SXW
              I4=SXX*SWW-SXW*SXW
              I5=SXZ*SXW-SXX*SZW
              I6=SXX*SZZ-SXZ*SXZ
              B1=0.
              B2=0.
              B3=0.
              IF (DET.GT.0.) B1=(I1*SXY+I2*SZY+I3*SWY)/DET
              IF (DET.GT.0.) B2=(I2*SXY+I4*SZY+I5*SWY)/DET
              IF (DET.GT.0.) B3=(I3*SXY+I5*SZY+I6*SWY)/DET
              SRES=0.
              SR1=0.
              SR2=0.
              SR3=0.
              SR4=0.
              SR5=0.
              SR6=0.
              DO 200 K=1,N
                   RES=Y(K)-B1*X(K)-B2*Z(K)-B3*W(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*X(K)*W(K)
                   SR4=SR4+RES*RES*Z(K)*Z(K)
                   SR5=SR5+RES*RES*Z(K)*W(K)
                   SR6=SR6+RES*RES*W(K)*W(K)
200          CONTINUE
             FN=FLOAT(N)
             VAR=SRES/(FN-3.)
             E1=0.
             E2=0.
             E3=0.
             E4=0.
             E5=0.
             E6=0.
             S1=0.
             S2=0.
             S3=0.
             R1=0.
             R2=0.
             R3=0.
             IF (DET.GT.0.) E1=FN*I1/DET
             IF (DET.GT.0.) E2=FN*I2/DET
             IF (DET.GT.0.) E3=FN*I3/DET
             IF (DET.GT.0.) E4=FN*I4/DET
             IF (DET.GT.0.) E5=FN*I5/DET
             IF (DET.GT.0.) E6=FN*I6/DET
             ES1=SR1/(FN-3.)
             ES2=SR2/(FN-3.)
             ES3=SR3/(FN-3.)
             ES4=SR4/(FN-3.)
             ES5=SR5/(FN-3.)
             ES6=SR6/(FN-3.)
             IF (DET.GT.0.) S1=SQRT(VAR*E1/FN)
             IF (DET.GT.0.) S2=SQRT(VAR*E4/FN)
             IF (DET.GT.0.) S3=SQRT(VAR*E6/FN)
           M1=E1*E1*ES1+E2*E2*ES4+E3*E3*ES6
           M2=E2*E1*ES2+E3*E1*ES3+E3*E2*ES5
           M3=E2*E2*ES1+E4*E4*ES4+E5*E5*ES6
           M4=E2*E4*ES2+E2*E5*ES3+E4*E5*ES5
           M5=E3*E3*ES1+E5*E5*ES4+E6*E6*ES6
           M6=E3*E5*ES2+E3*E6*ES3+E5*E6*ES5
          IF (DET.GT.0.) RR1=M1+2.*M2
          IF (DET.GT.0.) RR2=M3+2.*M4
          IF (DET.GT.0.) RR3=M5+2.*M6
          IF (DET.GT.0.) R1=SQRT(RR1/FN)
          IF (DET.GT.0.) R2=SQRT(RR2/FN)
          IF (DET.GT.0.) R3=SQRT(RR3/FN)
              RETURN
              END


       subroutine a(x,z,w,zz,ww,t1,t2,t3,t4,i,xg,zg,wg,b,m,kn,nobs)
               implicit none
               DOUBLE PRECISIon X(1),R,C,z(1),w(1),xd,wd,zd,zz(1),ww(1)
               double precision t1(1),t2(1),t3(1),t4(1)
               DOUBLE PRECISION W1,XG,zg,wg,r2,r3,r4
               integer j,nnc,nnd,kn,nobs,i
               double precision fnnc,fnnd,fkn,rro,f,f1,f2
               double precision dis(23972),ro(23972)
               double precision s1,s2,s3,s4,b,m,zzd,wwd
               R=T1(I)
               r2=t2(i)
               r3=t3(i)
               r4=t4(i)
               s1=15.1074
               s2=1.7918
               s3=1.2159876
               s4=1.7333
               C=0.
               xD=0.
               zd=0.
               wd=0.
               zzd=0.
               wwd=0.
               XG=0.
               zg=0.
               wg=0.
               b=0.
               m=0.
                   do 301 j=1,nobs
        f1=((t1(j)-r)/s1)**2+((t2(j)-r2)/s2)**2
        f2=((t3(j)-r3)/s3)**2+((t4(j)-r4)/s4)**2
        f=f1+f2
             dis(j)=sqrt(f)
                     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
                   if(j.eq.i) w1=0.
               endif
               xD=xD+W1*X(j)
               zD=zD+W1*z(j)
               wD=wD+W1*w(j)
               zzd=zzd+w1*zz(j)
               wwd=wwd+w1*ww(j)
303            CONTINUE
                XG=xD
                zG=zD
                wG=wD
                b=zzd
                m=wwd

              RETURN
               END


               subroutine aknn(x,t1,t2,t3,t4,i,xg,kn,nobs)
               implicit none
               DOUBLE PRECISIon X(1),R,C,xd
               double precision t1(1),t2(1),t3(1),t4(1)
               DOUBLE PRECISION W1,XG,r2,r3,r4
               integer j,nnc,nnd,kn,nobs,i
               double precision fnnc,fnnd,fkn,rro,f,f1,f2
               double precision dis(23972),ro(23972)
               double precision s1,s2,s3,s4
               R=T1(I)
               r2=t2(i)
               r3=t3(i)
               r4=t4(i)
               s1=15.1074
               s2=1.7918
               s3=1.2159876
               s4=1.7333
               C=0.
               xD=0.
               XG=0.
                   do 301 j=1,nobs
             f1=((t1(j)-r)/s1)**2+((t2(j)-r2)/s2)**2
             f2=((t3(j)-r3)/s3)**2+((t4(j)-r4)/s4)**2
             f=f1+f2
             dis(j)=sqrt(f)
                     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
                   if(j.eq.i) w1=0.
               endif
               xD=xD+W1*X(j)
303            CONTINUE
                XG=xD

              RETURN
               END

            SUBROUTINE SORT(N,RA)
            implicit none
            double precision ra(23972),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
