
*   /*    ----- 労働供給の4者択一モデル ----   */
*   /* Maximum Likelihood Estimation by Quadratic Hill-Climbing  */
*   /*  OUTPUT='LS4EST.TXT'                                      */
*   /*  OPTIONS:                                                 */
*   /*    1) LEVEL OF REPORTING                                  */
*   /*            LVRPU=(0,1,2)                                  */
*   /*              DEFAULT IS 0.                                */
*   /*    2) METHOD OF CALCULATING VARIANCE OF PARAMETERS        */
*   /*            BHHHVR=('0','1')                               */
*   /*                    '0': HESSIAN MATRIX                    */
*   /*                    '1': BHHH (Outer Product of Gradient)  */
*   /*              DEFAULT IS '1'.                              */
*   /*    3) METHOD OF CALCULATING HESSIAN MATRIX IN             */
*   /*         THE QUADRATIC HILL-CLIMBING                       */
*   /*            BHHHQH=('0','1')                               */
*   /*                    '0': DERIVATIVES IN SECOND ORDER       */
*   /*                    '1': BHHH (Outer Product of Gradient)  */
*   /*              DEFAULT IS '0'.                              */
*   /*    4) NUMBER OF TRIALS OF PARAMETER SEARCH                */
*   /*            NTRY=(1,2,...)                                 */
*   /*              DEFAULT IS 2.                                */
*   /*    5) MAXIMUM ITERATION                                   */
*   /*            MAXIT=(1,2,....)                               */
*   /*              DEFAULT IS 50.                               */
*   /*    6) ESTIMATION PERIODS STARTS AT NEN1 AND ENDS AT NEN2 .*/
*   /*            NEN1,NEN2=(1,2,..,6), NEN1<NEN2                */
*   /*              DEFAULT IS NEN1=1, NEN2=6                    */
*   /* "ITERATIVE" ESTIMATION OF PARAMETERS                      */
*   /*   BY QUADRATIC HILL-CLIMBING METHOD INCLUDING SQUEEZING.  */
*   /* ITERATION TERMINATES WHEN RATE>RATEMX OR RATES<RATEMN.    */

      IMPLICIT REAL*8(A-Z)
      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 MUOBS(L,4)
      REAL*8 MUEST(L,4)
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)
      INTEGER NPRM,NPRMS,LPRMS(M)
      INTEGER NSMPL
      INTEGER MAXIT,LVRPU
      REAL*8 PRMB(M)
      REAL*8 SIGMAV(M,M)
      INTEGER CODE
      INTEGER I,J,K,NSTEP,IT,NTRY,KOUNT,KTRY
      REAL*8 EPS,PI

      CHARACTER*1 BHHHVR
      CHARACTER*1 BHHHQH

      REAL*8 ESTPRM
      REAL*8 VARPRM
      REAL*8 SDPRM
      REAL*8 TVALUE
      CHARACTER*21 LBLPRM(M)
      CHARACTER*11 LBLSTA(M)


      COMMON/BLKOBN/OBSN(L,4)
      COMMON/BLKOB/MUOBS
      COMMON/BLKES/MUEST
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKOP/MAXIT,LVRPU
      COMMON/BLKMT/NSMPL
      COMMON/BLKPR/PRMB
      COMMON/BLKEPS/EPS
      COMMON/BLKPI/PI
      COMMON/BLKSW/BHHHQH

      DATA LBLPRM/'$\gamma_2$           '
     :           ,'$\gamma_3$           '
     :           ,'$\overline{\gamma_4}$'
     :           ,'$\gamma_5$           '
     :           ,'$\gamma_4^0$         '
     :           ,'$\sigma$             ',14*' '/
      DATA LBLSTA/'Estimate : '
     :           ,'Variance : '
     :           ,'S.D.     : '
     :           ,'T-value  : ',16*' '/


      EPS=0.1D-8
      PI =ACOS(-1.0D0)

      BHHHVR='1'
      BHHHQH='0'
      NTRY=2

      OPEN(UNIT=30,FILE='LS4EST.TXT')

      NSTEP=36
      KTRY=NTRY*NSTEP
      KOUNT=0
      CALL PREP4(NSTEP,CODE)
      IF(CODE.NE.0) STOP
      DO 30 IT=1,NTRY
       DO 20 K=1,NSTEP
        KOUNT=KOUNT+1
        CALL SETDRV(K,CODE)
        IF(CODE.NE.0) STOP
        CALL QHILML
        IF(KOUNT.EQ.KTRY)THEN
          IF( BHHHVR.EQ.'1' )THEN
            CALL VRBHHH(SIGMAV)
          ELSE
            CALL VAR(SIGMAV)
          END IF
          WRITE(6,'(1H ,A)')'----<< RESULT OF ESTIMATION >>----'
          WRITE(30,'(1H ,A)')'----<< RESULT OF ESTIMATION >>----'
          DO 110 I=1,NPRMS
            ESTPRM=PRMB(LPRMS(I))
            VARPRM=SIGMAV(I,I)
            SDPRM =SQRT(ABS(VARPRM))
            TVALUE=PRMB(LPRMS(I))/SDPRM
            WRITE(6,*)LBLPRM(LPRMS(I))
            WRITE(6,3100)LBLSTA(1),ESTPRM
            WRITE(6,3100)LBLSTA(2),VARPRM
            WRITE(6,3100)LBLSTA(3),SDPRM
            WRITE(6,3100)LBLSTA(4),TVALUE
            WRITE(30,*)LBLPRM(LPRMS(I))
            WRITE(30,3100)LBLSTA(1),ESTPRM
            WRITE(30,3100)LBLSTA(2),VARPRM
            WRITE(30,3100)LBLSTA(3),SDPRM
            WRITE(30,3100)LBLSTA(4),TVALUE
  110     CONTINUE
          WRITE(6,'(1H ,A)')
     :    '--<< OBSERVED PROB. VS. PREDICTED PROB. OF LABOR SUPPLY >>--'
          WRITE(30,'(1H ,A)')
     :    '--<< OBSERVED PROB. VS. PREDICTED PROB. OF LABOR SUPPLY >>--'
          DO 10 I=1,NSMPL
            WRITE(6,'(1H ,A,I5)')'----- GROUP:',I
            WRITE(6,'(1H ,A,4F10.6)')'MUOBS:',(MUOBS(I,J),J=1,4)
            WRITE(6,'(1H ,A,4F10.6)')'MUEST:',(MUEST(I,J),J=1,4)
            WRITE(30,'(1H ,A,I5)')'----- GROUP:',I
            WRITE(30,'(1H ,A,4F10.6)')'MUOBS:',(MUOBS(I,J),J=1,4)
            WRITE(30,'(1H ,A,4F10.6)')'MUEST:',(MUEST(I,J),J=1,4)
   10     CONTINUE
        ENDIF
   20  CONTINUE
   30 CONTINUE

      CLOSE(UNIT=30)

 3100 FORMAT(1H ,A,G20.13)

      STOP
      END
      BLOCKDATA
      INTEGER M
      PARAMETER(M=20)
      INTEGER NPRMSL(36),JPRMSL(36,6)
      CHARACTER TPRMS(M)*8
      COMMON/BLKDR/NPRMSL,JPRMSL
      COMMON/BLKTP/TPRMS
      DATA NPRMSL/20*3,15*2,6/
      DATA JPRMSL
     :/1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,4
     :  ,1,1,1,1,1,2,2,2,2,3,3,3,4,4,5, 1,
     : 2,2,2,2,3,3,3,4,4,5,3,3,3,4,4,5,4,4,5,5
     :  ,2,3,4,5,6,3,4,5,6,4,5,6,5,6,6, 2,
     : 3,4,5,6,4,5,6,5,6,6,4,5,6,5,6,6,5,6,6,6,15*0, 3,
     : 35*0, 4, 35*0, 5, 35*0, 6/
      DATA TPRMS/'G2','G3','G4','G5','G4O','SIGMA',14*' '/
      END
      SUBROUTINE PREP4(NSTEP,CODE)
      IMPLICIT REAL*8(A-Z)
      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 MUOBS(L,4)
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)
      REAL*8 WA(10),HA(10),VA(10)
      INTEGER NPRM,NPRMS,LPRMS(M)
      INTEGER NSMPL
      INTEGER MAXIT,LVRPU
      REAL*8 PRMB(M)
      CHARACTER TPRMS(M)*8
      INTEGER CODE
      INTEGER I,K,N
      INTEGER NEN,LNEN(6)
      INTEGER NEN1,NEN2
      CHARACTER MUFILE*10,FPARM*20

      INTEGER NPRMSL(36),JPRMSL(36,6),NSTEP

      COMMON/BLKOBN/OBSN(L,4)
      COMMON/BLKOB/MUOBS
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKOP/MAXIT,LVRPU
      COMMON/BLKMT/NSMPL
      COMMON/BLKPR/PRMB
      COMMON/BLKDR/NPRMSL,JPRMSL
      COMMON/BLKTP/TPRMS
      SAVE

      DATA MUFILE/'TABMUA.S??'/
      DATA NEN/6/
      DATA LNEN/46,49,52,54,57,62/


*  NPRM:totoal number of parameters.

      NPRM=6


*  MAXIT:maximum number of iterations.
*  LVRPU:print control variable(0=minimum,1=maximum).

      MAXIT=50
      LVRPU= 0


      NEN1=1
      NEN2=6

      FPARM='LS4PARM.TXT'
      WRITE(6,'(1H ,2A)')'>> input PARAMETER set file is ',FPARM
      WRITE(30,'(1H ,2A)')'>> input PARAMETER set file is ',FPARM

      OPEN(UNIT=10,FILE=FPARM,STATUS='OLD')
      READ(10,*)PRMB(1)
      READ(10,*)PRMB(2)
      READ(10,*)PRMB(3)
      READ(10,*)PRMB(4)
      READ(10,*)PRMB(5)
      READ(10,*)PRMB(6)
      CLOSE(UNIT=10)

      OPEN(UNIT=10,FILE='A60WV.DAT',STATUS='OLD')
      DO 110 N=1,NEN
        READ(10,'(5X,F10.3)')VA(N)
  110 CONTINUE
      CLOSE(UNIT=10)
      OPEN(UNIT=10,FILE='HWCENSUS.DAT',STATUS='OLD')
      DO 120 N=1,NEN
        READ(10,'(5X,2F10.3)')HA(N),WA(N)
  120 CONTINUE
      CLOSE(UNIT=10)

      K=0
      DO 130 N=NEN1,NEN2
        WRITE(MUFILE(9:10),'(I2)')LNEN(N)
        OPEN(UNIT=10,FILE=MUFILE,STATUS='OLD')
   10   CONTINUE
          K=K+1
          READ(10,1000,END=20)WEITAB(K),XITAB(K),(MUOBS(K,I),I=1,4)
          DO 140 I=1,4
            OBSN(K,I)=WEITAB(K)*MUOBS(K,I)
  140     CONTINUE
          WWTAB(K)=WA(N)
          HWTAB(K)=HA(N)
          VWTAB(K)=VA(N)
          GOTO 10
   20   CONTINUE
        K=K-1
        CLOSE(UNIT=10)
  130 CONTINUE
      NSMPL=K

      IF(MAXIT.LE.0)MAXIT=50
      WRITE(6,'(4X,A,I4)')
     :  'NUMBER OF GROUP=',NSMPL,'PARAMETERS=',NPRM,
     :  'MAXIMUM ITERATION=',MAXIT
      WRITE(30,'(4X,A,I4)')
     :  'NUMBER OF GROUP=',NSMPL,'PARAMETERS=',NPRM,
     :  'MAXIMUM ITERATION=',MAXIT

      WRITE(6,'(1H ,A/1H ,A,2I4)')
     :'<< ESTIMATION OF PARAMETERS BY QUADRATIC HILL CLIMBING METHOD >>'
     :,'----- 労働供給の4者択一モデル ----',LNEN(NEN1),LNEN(NEN2)
      WRITE(30,'(1H ,A/1H ,A,2I4)')
     :'<< ESTIMATION OF PARAMETERS BY QUADRATIC HILL CLIMBING METHOD >>'
     :,'----- 労働供給の4者択一モデル ----',LNEN(NEN1),LNEN(NEN2)
      WRITE(6,'(1H ,A)') '____ ITERATIVE SEARCH CONDITION ____'
      WRITE(6,'(1H ,A)')'INITIAL PARAMETERS'
      WRITE(6,'(1H ,6G13.6)')(PRMB(I),I=1,NPRM)
      WRITE(6,'(1H ,A)')'DRIVING PARAMETERS'
      WRITE(30,'(1H ,A)') '____ ITERATIVE SEARCH CONDITION ____'
      WRITE(30,'(1H ,A)')'INITIAL PARAMETERS'
      WRITE(30,'(1H ,6G13.6)')(PRMB(I),I=1,NPRM)
      WRITE(30,'(1H ,A)')'DRIVING PARAMETERS'
      DO 70 K=1,NSTEP
        WRITE(6,'(1H ,A,I2.2,A,13A8)')'STEP',K,':PARAMETERS TO DRIVE:',
     :    (TPRMS(JPRMSL(K,I)),I=1,NPRMSL(K))
        WRITE(30,'(1H ,A,I2.2,A,13A8)')'STEP',K,':PARAMETERS TO DRIVE:',
     :    (TPRMS(JPRMSL(K,I)),I=1,NPRMSL(K))
   70 CONTINUE

      CODE=0
      IF(NSMPL.GT.L.OR.NPRM.GT.M)THEN
        WRITE(6,'(1H ,A)')'++++++++  EXPAND THE MATRIX SIZE.'
        WRITE(6,'(1H ,4(A,I5,5X))')
     :  'L=',L,'NSMPL=',NSMPL,'M=',M,'NPRM=',NPRM
        CODE=1
      ENDIF

      RETURN
 1000 FORMAT(3X,F7.0,F10.2,4F10.7)
 2000 FORMAT(1H ,5G13.6,F9.6)
      END
      SUBROUTINE SETDRV(IDRV,CODE)
      IMPLICIT REAL*8(A-Z)
      INTEGER M
      PARAMETER(M=20)

*   FUNCTIONS
      INTEGER PRECHK

      CHARACTER LRES(5)*70

      INTEGER KR
      PARAMETER(KR=70)
      INTEGER KT,IDTR(KR)

      INTEGER IDRV,CODE
      REAL*8 PRMB(M)
      INTEGER NPRM,NPRMS,LPRMS(M)
      INTEGER MAXIT,LVRPU
      INTEGER NPRMSL(36),JPRMSL(36,6)
      CHARACTER TPRMS(M)*8,CPRMS(M)*1
      INTEGER I,J

      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKOP/MAXIT,LVRPU
      COMMON/BLKDR/NPRMSL,JPRMSL
      COMMON/BLKPR/PRMB
      COMMON/BLKTP/TPRMS
      SAVE

      DATA LRES/
     :'1=A1<0,1-2=SLOPE RIGHT DOWN,3=MU OF INCOME,4-6=MU OF LEISURE,',
     :'7-18=CONVEX TO ORIGINE,19=H0>IMAX,20=H2<0,',
     :'21=φ関数の勾配正,22-23=0<φ[H(d)=0]=f[H(d)=0],',
     :'24=φ[H(d)=0]=f[H(d)=0]<h,25=q1の座標<h,26=0<q1の座標,',
     :'27=q4の座標<a,28=h<q4の座標,29=q1の計算の平方根,30-31=G4>0'/

      DO 50 I=1,NPRM
        CPRMS(I)='0'
   50 CONTINUE
      NPRMS=NPRMSL(IDRV)
      DO 60 I=1,NPRMS
        J=JPRMSL(IDRV,I)
        LPRMS(I)=J
        CPRMS(J)='1'
   60 CONTINUE
      WRITE(6,'(1H ,A,I5)')'INITIAL PARAMETERS OF SEARCH PHASE:',IDRV
      WRITE(6,'(1H ,6G13.6)')(PRMB(I),I=1,NPRM)
      WRITE(6,'(1H ,A)') 'CONTROL DATA OF DRIVING PARAMETERS'
      WRITE(6,'(6(2X,A1))') (CPRMS(I),I=1,NPRM)
      WRITE(6,'(1H ,A,13A8)') 'PARAMETERS TO DRIVE:',
     :        (TPRMS(LPRMS(I)),I=1,NPRMS)
      WRITE(30,'(1H ,A,I5)')'INITIAL PARAMETERS OF SEARCH PHASE:',IDRV
      WRITE(30,'(1H ,6G13.6)')(PRMB(I),I=1,NPRM)
      WRITE(30,'(1H ,A)') 'CONTROL DATA OF DRIVING PARAMETERS'
      WRITE(30,'(6(2X,A1))') (CPRMS(I),I=1,NPRM)
      WRITE(30,'(1H ,A,13A8)') 'PARAMETERS TO DRIVE:',
     :        (TPRMS(LPRMS(I)),I=1,NPRMS)

      CODE=PRECHK(PRMB,LVRPU,IDTR,KT)
      IF(CODE.EQ.1)THEN
        WRITE(6,'(1H ,A)')
     :   '###### INITIAL PARAMETER SET VIOLATES TH.RES. ######'
        WRITE(30,'(1H ,A)')
     :   '###### INITIAL PARAMETER SET VIOLATES TH.RES. ######'
        IF(KT.GT.0)THEN
          WRITE(6,'(1H ,A)')LRES
          WRITE(30,'(1H ,A)')LRES
          WRITE(6,'(1H ,A)')'  ### RESTRICTION ID:'
          WRITE(6,'(20(1X,I2))')(IDTR(I),I=1,KT)
          WRITE(30,'(1H ,A)')'  ### RESTRICTION ID:'
          WRITE(30,'(20(1X,I2))')(IDTR(I),I=1,KT)
        ENDIF
      ENDIF

      RETURN
      END
      SUBROUTINE VAR(SIGMAV)
      IMPLICIT REAL*8(A-Z)

      INTEGER M
      PARAMETER(M=20)
      INTEGER NSMPL,NPRM,NPRMS,LPRMS(M)

      REAL*8 PARM(M)
      REAL*8 ED2LDP(M,M)
      REAL*8 EVEC(M),SOL(M)
      REAL*8 SIGMAV(M,M)
      REAL*8 XINFMT(M,M)
      INTEGER JJ,KK

      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      COMMON/BLKPR/PARM
      SAVE


      CALL HSSIAN(PARM,ED2LDP)

      DO 120 KK=1,NPRMS
        DO 110 JJ=1,NPRMS
          XINFMT(JJ,KK)=-ED2LDP(JJ,KK)
  110   CONTINUE
  120 CONTINUE

      DO 210 KK=1,NPRMS

        DO 220 JJ=1,NPRMS
          EVEC(JJ)=0.0D0
  220   CONTINUE
        EVEC(KK)=1.0D0

        CALL SLVLEQ(XINFMT,EVEC,NPRMS,SOL)

        DO 230 JJ=1,NPRMS
          SIGMAV(JJ,KK)=SOL(JJ)/DBLE(NSMPL)
  230   CONTINUE

  210 CONTINUE

      RETURN
      END
      SUBROUTINE VRBHHH(SIGMAV)

      IMPLICIT REAL*8(A-Z)

      INTEGER M
      PARAMETER(M=20)

      REAL*8 SIGMAV(M,M)

      INTEGER LPRMS(M)
      INTEGER NSMPL,NPRM,NPRMS

      REAL*8 PARM(M)
      REAL*8 SUMOPG(M,M)
      REAL*8 EVEC(M)
      REAL*8 SOL(M)
      REAL*8 XINFMT(M,M)

      INTEGER JJ,KK

      REAL*8 SMOBSN


      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      COMMON/BLKPR/PARM
      SAVE


      CALL OPGSUM(PARM,SUMOPG,SMOBSN)

      DO 110 KK=1,NPRM
        DO 120 JJ=1,NPRM
          XINFMT(JJ,KK)=SUMOPG(JJ,KK)/SMOBSN
  120   CONTINUE
  110 CONTINUE

      DO 210 KK=1,NPRM

        DO 220 JJ=1,NPRM
          EVEC(JJ)=0.0D0
  220   CONTINUE
        EVEC(KK)=1.0D0

        CALL SLVLEQ(XINFMT,EVEC,NPRM,SOL)

        DO 230 JJ=1,NPRM
          SIGMAV(JJ,KK)=SOL(JJ)
  230   CONTINUE

  210 CONTINUE 

      RETURN
      END
      SUBROUTINE OPGSUM(PARM,SUMOPG,SMOBSN)

      IMPLICIT REAL*8(A-Z)

      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 PARM(M)
      REAL*8 SUMOPG(M,M)
      REAL*8 SMOBSN

      REAL*8 DLGLLP(L,4,M)  
      REAL*8 OPGLGL(L,4,M,M)

      REAL*8 OBSN (L,4)
      REAL*8 MUOBS(L,4)
      REAL*8 MUEST(L,4)
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)

      INTEGER NSMPL,NPRM
      INTEGER LPRMS(M)
      INTEGER NPRMS
      INTEGER I,J,IP,JP

      REAL*8 WW,HW,VW,XI
      REAL*8 DPDPRM(1:4,M)

      COMMON/BLKOBN/OBSN
      COMMON/BLKOB/MUOBS
      COMMON/BLKES/MUEST
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      SAVE


      SMOBSN=0.0D0

      DO 10 IP=1,NPRM
        DO 20 JP=1,NPRM
          SUMOPG(IP,JP)=0.0D0
   20   CONTINUE
   10 CONTINUE


      DO 110 I=1,NSMPL

        WW=WWTAB(I)
        HW=HWTAB(I)
        VW=VWTAB(I)
        XI=XITAB(I)

        CALL NABLAP(DPDPRM,PARM,WW,HW,VW,XI)

        DO 120 J=1,2
          DO 130 IP=1,NPRM
            DLGLLP(I,J,IP)=OBSN(I,J)*DPDPRM(J,IP)/MUEST(I,J)
            SMOBSN=SMOBSN+OBSN(I,J)
  130     CONTINUE
  120   CONTINUE


        DO 140 J=1,2

          DO 150 IP=1,NPRM
            DO 160 JP=1,NPRM
              OPGLGL(I,J,IP,JP)=DLGLLP(I,J,IP)
     :                         *DLGLLP(I,J,JP)
  160       CONTINUE
  150     CONTINUE

          DO 170 IP=1,NPRM
            DO 180 JP=1,NPRM
              SUMOPG(IP,JP)=SUMOPG(IP,JP)+OPGLGL(I,J,IP,JP)
  180       CONTINUE
  170     CONTINUE

  140   CONTINUE

  110 CONTINUE

      RETURN
      END
      SUBROUTINE QHILML
      IMPLICIT REAL*8(A-Z)

      INTEGER L,M
      PARAMETER(L=200,M=20)

      INTEGER KR
      PARAMETER(KR=70)
      INTEGER KT,IDTR(KR)
      CHARACTER LRES(5)*70

*   FUNCTIONS
      INTEGER PRECHK
      REAL*8 SUMLGL

      REAL*8 MUOBS(L,4)
      INTEGER NSMPL,NPRM
      INTEGER LVRPU
      INTEGER NPRMS,LPRMS(M)

      REAL*8 F(L,4),MTRXA(L,4,M)
      REAL*8 FC(L,4),FCM(L,4),FCP(L,4)
      INTEGER ISQ,ICSQ,CODE
      REAL*8 SQRATE
      CHARACTER EOSQ*1

      INTEGER IT,MAXIT,I,J,MAXSQ
      INTEGER NROT
      REAL*8 EPSFLT
      REAL*8 PARM(M)
      REAL*8 SFXHES(M,M)
      REAL*8 DSDPRM(M)
      REAL*8 PARM2(M),PARMP(M),PARMM(M)
      REAL*8 PDELTA(M)
      REAL*8 EVAL(M),EVEC(M,M)
      REAL*8 LMDMIN,EVECMN(M)
      REAL*8 ALPHA,NORMDF
      REAL*8 F0,F1,FP,FM
      REAL*8 RATE
      REAL*8 RATES
      REAL*8 RINIT
      REAL*8 MRATE
      CHARACTER*1 FLAG,EOS,RED

      REAL*8 RATEMX,RATEMN,RSHIFT,RMAXP
      CHARACTER*1 SW

      REAL*8 SUMOPG(M,M)
      REAL*8 SMOBSN

      CHARACTER*1 BHHHQH

      COMMON/BLKOBN/OBSN(L,4)
      COMMON/BLKOB/MUOBS
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKOP/MAXIT,LVRPU
      COMMON/BLKMT/NSMPL
      COMMON/BLKPR/PARM
      COMMON/BLKSW/BHHHQH
      SAVE

      DATA LRES/
     :'1=A1<0,1-2=SLOPE RIGHT DOWN,3=MU OF INCOME,4-6=MU OF LEISURE,',
     :'7-18=CONVEX TO ORIGINE,19=H0>IMAX,20=H2<0,',
     :'21=φ関数の勾配正,22-23=0<φ[H(d)=0]=f[H(d)=0],',
     :'24=φ[H(d)=0]=f[H(d)=0]<h,25=q1の座標<h,26=0<q1の座標,',
     :'27=q4の座標<a,28=h<q4の座標,29=q1の計算の平方根,30-31=G4>0'/

      EOS='0'
      EPSFLT=0.1D-12
      RATE=1.0D-6
      RATES=1.0D0
      RINIT=1.0D-6
      MRATE=SQRT(10.0D0)

      MAXSQ=20
      SQRATE=0.9D0

      RATEMX=0.1D20
      RATEMN=1.0D0/RATEMX
      RSHIFT=0.1D-8
      SW='0'

      DO 10 I=1,NPRM
        PARM2(I)=PARM(I)
   10 CONTINUE
      CALL LOGL(F,PARM)
      F0=SUMLGL(F,NSMPL)

      DO 100 IT=0,MAXIT

        WRITE(6,'(1H ,A,I5)')'-------- ITERATION=',IT
        WRITE(6,1000)(PARM(I),I=1,NPRM)
        IF(LVRPU.GT.0)THEN
          WRITE(30,'(1H ,A,I5)')'-------- ITERATION=',IT
          WRITE(30,1000)(PARM(I),I=1,NPRM)
        ENDIF

        CALL LOGL(F,PARM)
        CALL JACOB(PARM,MTRXA)
        CALL NABLAL(MTRXA,DSDPRM)
        IF( BHHHQH.EQ.'1' )THEN
          CALL OPGSUM(PARM,SUMOPG,SMOBSN)
          DO 20 I=1, NPRMS
            DO 30 J=1, NPRMS
              SFXHES(I,J)=-SUMOPG(LPRMS(I),LPRMS(J))
   30       CONTINUE
   20     CONTINUE
        ELSE
          CALL HSSIAN(PARM,SFXHES)
        ENDIF
        CALL JACOBI(SFXHES,NPRMS,EVAL,EVEC,NROT,FLAG)

        WRITE(6,'(1H ,A)')'NABLA LOG L'
        WRITE(6,'(1X,6G12.4)')(-DSDPRM(I),I=1,NPRMS)
        WRITE(6,'(1H ,A)')'EIGEN VALUES'
        WRITE(6,'(1X,6G12.4)')(-EVAL(I),I=1,NPRMS)
        WRITE(6,'(1H ,A,G20.13)')'LOG OF LILELIHOOD FUNCTION=',-F0

        IF(LVRPU.GT.0)THEN
          WRITE(30,'(1H ,A)')'NABLA LOG L'
          WRITE(30,'(1X,6G12.4)')(-DSDPRM(I),I=1,NPRMS)
          WRITE(30,'(1H ,A)')'EIGEN VALUES'
          WRITE(30,'(1X,6G12.4)')(-EVAL(I),I=1,NPRMS)
          WRITE(30,'(1H ,A,G20.13)')'LOG OF LIKELIHOOD FUNCTION=',-F0
        ENDIF

        NORMDF=0.0
        LMDMIN=+1.0D+20
        DO 200 I=1,NPRMS
          IF(EVAL(I).LT.LMDMIN)THEN
            LMDMIN=EVAL(I)
            DO 210 J=1,NPRMS
              EVECMN(J)=EVEC(J,I)
  210       CONTINUE
          END IF
          NORMDF=NORMDF+DSDPRM(I)**2
  200   CONTINUE
        NORMDF=SQRT(NORMDF)
        IF(NORMDF.GE.EPSFLT)THEN
          SW='0'
          RMAXP=RATE
  300     CONTINUE
            ALPHA=LMDMIN-RATE*NORMDF
            IF(ALPHA.LT.0.0D0)THEN
              DO 310 I=1,NPRMS
                SFXHES(I,I)=SFXHES(I,I)-ALPHA
  310         CONTINUE
            END IF

            CALL SLVLEQ(SFXHES,DSDPRM,NPRMS,PDELTA)

            EOSQ='0'
            DO 600 ISQ=0,MAXSQ
              ICSQ=ISQ
              IF(ISQ.GT.0)THEN
                DO 610 I=1,NPRMS
                  PDELTA(I)=PDELTA(I)*SQRATE**ISQ
  610           CONTINUE
              ENDIF
              DO 620 I=1,NPRMS
                PARM2(LPRMS(I))=PARM(LPRMS(I))-PDELTA(I)
  620         CONTINUE

              CODE=PRECHK(PARM2,LVRPU,IDTR,KT)
              IF(CODE.EQ.0)EOSQ='1'
              IF(EOSQ.EQ.'1')GOTO 630
  600       CONTINUE
            IF(EOSQ.EQ.'0')THEN
              WRITE(6,'(1H ,2A,I4)')
     :         '### VIOLATION OF THEORETICAL RESTRICTION.'
     :         ,'SEARCH WILL TERMINAME. MAXSQ=',MAXSQ
              WRITE(6,'(1H ,A)')' DELTA OF PARAMETERS:'
              WRITE(6,'(6(1X,G12.5))')(PDELTA(I),I=1,NPRMS)
              WRITE(30,'(1H ,2A,I4)')
     :         '### VIOLATION OF THEORETICAL RESTRICTION.'
     :         ,'SEARCH WILL TERMINAME. MAXSQ=',MAXSQ
              WRITE(30,'(1H ,A)')' DELTA OF PARAMETERS:'
              WRITE(30,'(6(1X,G12.5))')(PDELTA(I),I=1,NPRMS)
              IF(KT.GT.0)THEN
                WRITE(6,'(1H ,A)')'  ### RESTRICTION ID:'
                WRITE(6,'(20(1X,I2))')(IDTR(I),I=1,KT)
                WRITE(30,'(1H ,A)')'  ### RESTRICTION ID (WIF.):'
                WRITE(30,'(20(1X,I2))')(IDTR(I),I=1,KT)
                WRITE(6,'(1H ,A)')LRES
                WRITE(30,'(1H ,A)')LRES
              ENDIF
              GOTO 500
            ENDIF

  630       CONTINUE
            IF(EOSQ.EQ.'1')THEN
              CALL LOGL(FC,PARM2)
              F1=SUMLGL(FC,NSMPL)
            ENDIF

            RED='0'
            WRITE(6,'(1H ,A,G13.6)')'     DF^=0: R=',RATE
            WRITE(6,'(1H ,A,G20.13)')'        LOG L=',-F1
            IF(LVRPU.GT.0)THEN
              WRITE(30,'(1H ,A,G13.6)')'     DF^=0: R=',RATE
              WRITE(30,'(1H ,A,G20.13)')'        LOG L=',-F1
            ENDIF
            IF(F1.LT.F0)THEN
              RED='1'
              DO 350 I=1,NPRMS
                PARM(LPRMS(I))=PARM2(LPRMS(I))
  350         CONTINUE
              DO 370 I=1,NSMPL
                DO 360 J=1,4
                  F(I,J)=FC(I,J)
  360           CONTINUE
  370         CONTINUE
              F0=F1
              RATE=RATE/MRATE
            ELSE
              RATE=RATE*MRATE
            END IF
            IF(SW.EQ.'0')THEN
              IF(ABS(RATE).GT.RATEMX)THEN
                SW='1'
                RATE=RINIT*RSHIFT
                WRITE(6,'(1H ,A,D13.6)')'RESTART WITH R=',RATE
                WRITE(30,'(1H ,A,D13.6)')'RESTART WITH R=',RATE
              ELSE
              ENDIF
            ELSE
              IF(ABS(RATE).GT.RMAXP)THEN
                WRITE(6,'(1H ,A)')'SEARCH RESTART FAILED.'
                WRITE(30,'(1H ,A)')'SEARCH RESTART FAILED.'
                GOTO 500
              ELSE
                RATE=RATE*MRATE
              ENDIF
            ENDIF
          IF(.NOT.(EOSQ.EQ.'1'.AND.RED.EQ.'1'))GOTO 300
        ELSE
          IF(LMDMIN.LT.0.0D0)THEN
            WRITE(6,'(1H ,A)')'  ** SADDLE POINT **'
            WRITE(6,'(1H ,A)')'      GRADIENT:'
            WRITE(6,'(1X,6G12.4)')(DSDPRM(I),I=1,NPRMS)
            WRITE(6,'(1H ,A,G15.6)')'  MINIMUM EIGEN VALUE:',LMDMIN
            WRITE(30,'(1H ,A)')'  ** SADDLE POINT **'
            WRITE(30,'(1H ,A)')'      GRADIENT:'
            WRITE(30,'(1X,6G12.4)')(DSDPRM(I),I=1,NPRMS)
            WRITE(30,'(1H ,A,G15.6)')'  MINIMUM EIGEN VALUE:',LMDMIN
  400       CONTINUE
              DO 410 I=1,NPRM
                PARMM(I)=PARM(I)
                PARMP(I)=PARM(I)
  410         CONTINUE
              DO 420 I=1,NPRMS
                PARMM(LPRMS(I))=PARM(LPRMS(I))-RATES*EVECMN(I)
                PARMP(LPRMS(I))=PARM(LPRMS(I))+RATES*EVECMN(I)
  420         CONTINUE
              EOSQ='0'
              IF(PRECHK(PARMM,LVRPU,IDTR,KT).EQ.0.AND.
     :           PRECHK(PARMP,LVRPU,IDTR,KT).EQ.0)THEN
                RED='0'
                EOSQ='1'
                CALL LOGL(FCM,PARMM)
                FM=SUMLGL(FCM,NSMPL)
                CALL LOGL(FCP,PARMP)
                FP=SUMLGL(FCP,NSMPL)
                WRITE(6,'(1H ,A,G13.6)')'        R=',RATES
                WRITE(6,'(1H ,A,2G15.6)')'    LOG L=',-FM,-FP
                IF(LVRPU.GT.0)THEN
                  WRITE(30,'(1H ,A,G13.6)')'        R=',RATES
                  WRITE(30,'(1H ,A,2G15.6)')'    LOG L=',-FM,-FP
                ENDIF
                IF(FP.LT.F0)THEN
                  RED='1'
                  DO 430 I=1,NPRMS
                    PARM(LPRMS(I))=PARMP(LPRMS(I))
  430             CONTINUE
                  DO 450 I=1,NSMPL
                    DO 440 J=1,4
                      F(I,J)=FCP(I,J)
  440               CONTINUE
  450             CONTINUE
                  F0=FP
                ELSE
                  IF(FM.LT.F0)THEN
                    RED='1'
                    DO 460 I=1,NPRMS
                      PARM(LPRMS(I))=PARMM(LPRMS(I))
  460               CONTINUE
                    DO 480 I=1,NSMPL
                      DO 470 J=1,4
                        F(I,J)=FCM(I,J)
  470                 CONTINUE
  480               CONTINUE
                    F0=FM
                  ELSE
                    RATES=RATES/MRATE
                  END IF
                END IF
              ELSE
                RATES=RATES/MRATE
              END IF
            IF(.NOT.((EOSQ.EQ.'1'.AND.RED.EQ.'1')
     :           .OR.ABS(RATES).LT.RATEMN))GOTO 400
            IF(EOSQ.EQ.'0')THEN
              WRITE(6,'(1H ,A)')'SQUEEZING FAILED.'
              WRITE(30,'(1H ,A)')'SQUEEZING FAILED.'
              GOTO 500
            ELSE
              RATE=RINIT
              RATES=1.0D0
            ENDIF
          ELSE
            EOS='1'
          END IF
        END IF
        IF(EOS.EQ.'1')GOTO 500
  100 CONTINUE
      WRITE(6,'(1H ,A)')'$$$$$$ ITERATION EXAUSTED. MAXIT=',MAXIT
      WRITE(30,'(1H ,A)')'$$$$$$ ITERATION EXAUSTED. MAXIT=',MAXIT

  500 CONTINUE
      IF(EOS.EQ.'1')WRITE(6,'(1H ,A)')'CONVERGED.'
      WRITE(6,'(1H ,A)')'     GRADIENT:'
      WRITE(6,'(1X,6G12.4)')(-DSDPRM(I),I=1,NPRMS)
      WRITE(6,'(1H ,A)')'     EIGEN VALUES:'
      WRITE(6,'(1X,6G12.4)')(-EVAL(I),I=1,NPRMS)
      WRITE(6,'(1H ,A)')' ---- PARAMETERS:'
      WRITE(6,2000)(PARM(I),I=1,NPRM)
      WRITE(6,'(1H ,A,G20.13)')'   LOG OF LIKELIHOOD FUNCTION=',-F0
      IF(EOS.EQ.'1')WRITE(30,'(1H ,A)')'CONVERGED.'
      WRITE(30,'(1H ,A)')'     GRADIENT:'
      WRITE(30,'(1X,6G12.4)')(-DSDPRM(I),I=1,NPRMS)
      WRITE(30,'(1H ,A)')'     EIGEN VALUES:'
      WRITE(30,'(1X,6G12.4)')(-EVAL(I),I=1,NPRMS)
      WRITE(30,'(1H ,A)')' ---- PARAMETERS:'
      WRITE(30,2000)(PARM(I),I=1,NPRM)
      WRITE(30,'(1H ,A,G20.13)')'   LOG OF LIKELIHOOD FUNCTION=',-F0

      RETURN
 1000 FORMAT(1H ,5G13.6,F9.6)
 2000 FORMAT(1H ,G20.13)
      END
      SUBROUTINE LOGL(FLOGL,PARM)
      IMPLICIT REAL*8(A-Z)

      EXTERNAL LS4

      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 MUOBS(L,4)
      REAL*8 MUEST(L,4)
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)
      INTEGER NSMPL,NPRM
      INTEGER NPRMS,LPRMS(M)
      REAL*8 PARM(M)
      REAL*8 FLOGL(L,4)
      INTEGER LMU,II
      REAL*8 SGMW
      REAL*8 GW(6)
      REAL*8 WW,HW,VW,XI
      REAL*8 MUHAT(4)
      INTEGER ICON
      REAL*8 U0,U1,U4

      COMMON/BLKOBN/OBSN(L,4)
      COMMON/BLKOB/MUOBS
      COMMON/BLKES/MUEST
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      SAVE

      GW(1)=-1.0D0
      GW(2)=PARM(1)
      GW(3)=PARM(2)
      GW(4)=PARM(3)
      GW(5)=PARM(4)
      GW(6)=PARM(5)
      SGMW =PARM(6)

      DO 100 II=1,NSMPL
        WW=WWTAB(II)
        HW=HWTAB(II)
        VW=VWTAB(II)
        XI=XITAB(II)

        CALL LS4(MUHAT,GW,SGMW,WW,HW,VW,XI,ICON,U0,U1,U4)
        IF(ICON.NE.0)THEN
           WRITE(6,*)'ICON=',ICON
           RETURN
        ENDIF

        DO 200 LMU=1,4
          IF(ABS(MUHAT(LMU)).LT.0.1D-5)MUHAT(LMU)=0.1D-5
          FLOGL(II,LMU)=OBSN(II,LMU)*LOG(MUHAT(LMU))
          MUEST(II,LMU)=MUHAT(LMU)
  200   CONTINUE
  100 CONTINUE
      RETURN
      END
      FUNCTION SUMLGL(F,LV)
      IMPLICIT REAL*8(A-Z)
      INTEGER L
      PARAMETER(L=200)
      REAL*8 F(L,4),SUMLGL,SUM
      INTEGER LV,I,J

      SUM=0.0D0
      DO 120 I=1,LV
        DO 110 J=1,4
          SUM=SUM-F(I,J)
  110   CONTINUE
  120 CONTINUE
      SUMLGL=SUM

      RETURN
      END
      SUBROUTINE JACOB(PARM,MTRXA)
      IMPLICIT REAL*8(A-Z)

      EXTERNAL NABLAP

      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 PARM(M)
      REAL*8 MTRXA(L,4,M)
      REAL*8 MUOBS(L,4)
      REAL*8 MUEST(L,4)
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)
      INTEGER NSMPL,NPRM
      INTEGER NPRMS,LPRMS(M)
      INTEGER I,J,LX
      REAL*8 WW,HW,VW,XI
      REAL*8 DPDPRM(1:4,M)

      COMMON/BLKOBN/OBSN(L,4)
      COMMON/BLKOB/MUOBS
      COMMON/BLKES/MUEST
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      SAVE

      DO 200 I=1,NSMPL
        WW=WWTAB(I)
        HW=HWTAB(I)
        VW=VWTAB(I)
        XI=XITAB(I)

        CALL NABLAP(DPDPRM,PARM,WW,HW,VW,XI)

        DO 400 J=1,4
          DO 500 LX=1,NPRMS
            MTRXA(I,J,LX)=OBSN(I,J)*DPDPRM(J,LPRMS(LX))/MUEST(I,J)
  500     CONTINUE
  400   CONTINUE
  200 CONTINUE

      RETURN
      END
      SUBROUTINE NABLAL(MTRXA,DSDPRM)
      IMPLICIT REAL*8(A-Z)
      INTEGER L,M
      PARAMETER(L=200,M=20)

      REAL*8 DSDPRM(M)
      REAL*8 MTRXA(L,4,M)
      INTEGER NSMPL,NPRM
      INTEGER NPRMS,LPRMS(M)
      INTEGER I,J,K

      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      SAVE

      DO 110 K=1,NPRMS
        DSDPRM(K)=0.0D0
  110 CONTINUE

      DO 230 I=1,NSMPL
        DO 220 J=1,4
          DO 210 K=1,NPRMS
            DSDPRM(K)=DSDPRM(K)-MTRXA(I,J,K)
  210     CONTINUE
  220   CONTINUE
  230 CONTINUE

      RETURN
      END
      SUBROUTINE HSSIAN(PARM,SFXHES)
      IMPLICIT REAL*8(A-Z)

      INTEGER L,M,I,J
      PARAMETER(L=200,M=20)

      REAL*8   PARM(M)
      REAL*8   SFXHES(M,M)
      INTEGER  NPRM,NPRMS,LPRMS(M)
      REAL*8   PARMM(M),FM(L,4),MTRXM(L,4,M),DSDPM(M)
      REAL*8   PARMP(M),FP(L,4),MTRXP(L,4,M),DSDPP(M)
      REAL*8   HABA

      INTEGER MAXIT
      INTEGER LVRPU
      INTEGER KR
      PARAMETER(KR=70)
      INTEGER IDTR(KR)
      INTEGER KT
      INTEGER MAXSQZ
      INTEGER ISQZ
      INTEGER FLGSQZ

      REAL*8 DLTINI
      REAL*8 DLTRAT

      EXTERNAL PRECHK
      INTEGER PRECHK

      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKOP/MAXIT,LVRPU
      SAVE


      DLTINI = 0.1D-2
      DLTRAT = 0.8D0
      MAXSQZ = 100

      FLGSQZ = 0
      HABA = DLTINI

      DO 210 ISQZ = 1, MAXSQZ
        DO 220 I=1,NPRMS
          DO 230 J=1,NPRM
            PARMM(J)=PARM(J)
            PARMP(J)=PARM(J)
  230     CONTINUE
          PARMM(LPRMS(I))=PARM(LPRMS(I))*(1.0D0-HABA)-HABA
          PARMP(LPRMS(I))=PARM(LPRMS(I))*(1.0D0+HABA)+HABA
          IF(PRECHK(PARMM,LVRPU,IDTR,KT).EQ.0 .AND.
     :       PRECHK(PARMP,LVRPU,IDTR,KT).EQ.0)THEN
            FLGSQZ = 1
            GOTO 300
          ELSE
            HABA = HABA * DLTRAT
          END IF
  220   CONTINUE
  210 CONTINUE

      IF(FLGSQZ.EQ.0)THEN
        PRINT *, '  HESSIAN: SQUEEZING HABA FAILED.'
        WRITE(30,*) '  HESSIAN: SQUEEZING HABA FAILED.'
        STOP
      END IF

  300 CONTINUE

      DO 100 I=1,NPRMS
        DO 110 J=1,NPRM
          PARMM(J)=PARM(J)
          PARMP(J)=PARM(J)
  110   CONTINUE
        PARMM(LPRMS(I))=PARM(LPRMS(I))*(1.0D0-HABA)-HABA
        PARMP(LPRMS(I))=PARM(LPRMS(I))*(1.0D0+HABA)+HABA
        CALL LOGL(FM,PARMM)
        CALL LOGL(FP,PARMP)
        CALL JACOB(PARMM,MTRXM)
        CALL JACOB(PARMP,MTRXP)
        CALL NABLAL(MTRXM,DSDPM)
        CALL NABLAL(MTRXP,DSDPP)
        DO 120 J=1,NPRMS
          SFXHES(J,I)=(DSDPP(J)-DSDPM(J))
     :             /(2.0D0*PARM(LPRMS(I))*HABA+2.0D0*HABA)
  120   CONTINUE
  100   CONTINUE
      RETURN
      END
      SUBROUTINE SLVLEQ(A,B,N,X)
      IMPLICIT REAL*8(A-H,O-Z)
*   /*                                                        */
*   /*    SLVLEQ SOLVES LINEAR EQUATIONS BY SINGULAR          */
*   /*                                      VALUE             */
*   /*                                      DECOMPOSITION.    */
*   /*    REAL*8 VERSION                                      */
*   /*    CAUTION:IF W(*)<WMAX*X  THEN W(*) IS SET TO 0,      */
*   /*               AND X MAYBE AFFECT THE SOLUTION.         */
*   /*               X IS SET TO 1.0E-12,                     */
*   /*               AND IT MAY BE BETTER ABS(X) IS SMALL.    */
*   /*                                                        */
      INTEGER MP,NP
      PARAMETER(MP=20,NP=MP)

      INTEGER  N,J
      REAL*8 A(MP,NP),U(MP,NP),W(NP),V(NP,NP),B(NP),X(NP),WMAX,WMIN

      CALL SVDCMP(A,U,N,N,W,V)
      WMAX=0.0D0
      DO 130 J=1,N
        IF(W(J).GT.WMAX) WMAX=W(J)
  130 CONTINUE
      WMIN=WMAX*1.0D-12
      DO 140 J=1,N
        IF(W(J).LT.WMIN) W(J)=0.0D0
  140 CONTINUE
      CALL SVBKSB(U,W,V,N,N,B,X)
      RETURN
      END
      SUBROUTINE SVDCMP(A,U,M,N,W,V)
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER MP,NP
      PARAMETER(MP=20,NP=MP)

      INTEGER I,J,K,L,M,N,ITS,NM,JJ,LL
      REAL*8   A(MP,NP),U(MP,NP),W(NP),V(NP,NP),RV1(100)
      REAL*8   G,S,C,SCALE,H,F,ANORM,X,Y,Z

      IF(M.LT.N)THEN
        WRITE(6,'(1H ,A)')'YOU MUST AUGMENT A WITH EXTRA ZERO ROWS.'
        RETURN
      ENDIF

      DO 10 I=1,M
        DO 20 J=1,N
          U(I,J)=A(I,J)
   20   CONTINUE
   10 CONTINUE

      G=0.0D0
      SCALE=0.0D0
      ANORM=0.0D0
      DO 250 I=1,N
        L=I+1
        RV1(I)=SCALE*G
        G=0.0D0
        S=0.0D0
        SCALE=0.0D0
        IF(I.LE.M)THEN
          DO 110 K=I,M
            SCALE=SCALE+ABS(U(K,I))
  110     CONTINUE
          IF(SCALE.NE.0.0D0)THEN
            DO 120 K=I,M
              U(K,I)=U(K,I)/SCALE
              S=S+U(K,I)*U(K,I)
  120       CONTINUE
            F=U(I,I)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            U(I,I)=F-G
            IF(I.NE.N)THEN
              DO 150 J=L,N
                S=0.0D0
                DO 130 K=I,M
                  S=S+U(K,I)*U(K,J)
  130           CONTINUE
                F=S/H
                DO 140 K=I,M
                  U(K,J)=U(K,J)+F*U(K,I)
  140           CONTINUE
  150         CONTINUE
            ENDIF
            DO 160 K=I,M
              U(K,I)=SCALE*U(K,I)
  160       CONTINUE
          ENDIF
        ENDIF
        W(I)=SCALE*G
        G=0.0D0
        S=0.0D0
        SCALE=0.0D0
        IF((I.LE.M).AND.(I.NE.N))THEN
          DO 170 K=L,N
            SCALE=SCALE+ABS(U(I,K))
  170     CONTINUE
          IF(SCALE.NE.0.0D0)THEN
            DO 180 K=L,N
              U(I,K)=U(I,K)/SCALE
              S=S+U(I,K)*U(I,K)
  180       CONTINUE
            F=U(I,L)
            G=-SIGN(SQRT(S),F)
            H=F*G-S
            U(I,L)=F-G
            DO 190 K=L,N
              RV1(K)=U(I,K)/H
  190       CONTINUE
            IF(I.NE.M)THEN
              DO 230 J=L,M
              S=0.0D0
                DO 210 K=L,N
                  S=S+U(J,K)*U(I,K)
  210           CONTINUE
                DO 220 K=L,N
                  U(J,K)=U(J,K)+S*RV1(K)
  220           CONTINUE
  230         CONTINUE
            ENDIF
            DO 240 K=L,N
              U(I,K)=SCALE*U(I,K)
  240       CONTINUE
          ENDIF
        ENDIF
        ANORM=MAX(ANORM,(ABS(W(I))+ABS(RV1(I))))
  250 CONTINUE

      DO 320 I=N,1,-1
        IF(I.LT.N)THEN
          IF(G.NE.0.0D0)THEN
            DO 260 J=L,N
              V(J,I)=(U(I,J)/U(I,L))/G
  260       CONTINUE
            DO 290 J=L,N
              S=0.0D0
              DO 270 K=L,N
                S=S+U(I,K)*V(K,J)
  270         CONTINUE
              DO 280 K=L,N
                V(K,J)=V(K,J)+S*V(K,I)
  280         CONTINUE
  290       CONTINUE
          ENDIF
          DO 310 J=L,N
            V(I,J)=0.0D0
            V(J,I)=0.0D0
  310     CONTINUE
        ENDIF
        V(I,I)=1.0D0
        G=RV1(I)
        L=I
  320 CONTINUE

      DO 390 I=N,1,-1
        L=I+1
        G=W(I)
        IF(I.LT.N)THEN
          DO 330 J=L,N
            U(I,J)=0.0D0
  330     CONTINUE
        ENDIF
        IF(G.NE.0.0D0)THEN
          G=1.0D0/G
          IF(I.NE.N)THEN
            DO 360 J=L,N
              S=0.0D0
              DO 340 K=L,M
                S=S+U(K,I)*U(K,J)
  340         CONTINUE
              F=(S/U(I,I))*G
              DO 350 K=I,M
                U(K,J)=U(K,J)+F*U(K,I)
  350         CONTINUE
  360       CONTINUE
          ENDIF
          DO 370 J=I,M
            U(J,I)=U(J,I)*G
  370     CONTINUE
        ELSE
          DO 380 J=I,M
            U(J,I)=0.0D0
  380     CONTINUE
        ENDIF
        U(I,I)=U(I,I)+1.0D0
  390 CONTINUE

      DO 490 K=N,1,-1
        DO 480 ITS=1,30
          DO 410 LL=K,1,-1
            L=LL
            NM=L-1
            IF((ABS(RV1(L))+ANORM).EQ.ANORM) GOTO 2000
            IF((ABS(W(NM))+ANORM).EQ.ANORM) GOTO 1000
  410     CONTINUE
 1000     CONTINUE
          C=0.0D0
          S=1.0D0
          DO 430 I=L,K
            F=S*RV1(I)
            IF((ABS(F)+ANORM).NE.ANORM)THEN
              G=W(I)
              H=SQRT(F*F+G*G)
              W(I)=H
              H=1.0D0/H
              C= (G*H)
              S=-(F*H)
              DO 420 J=1,M
                Y=U(J,NM)
                Z=U(J,I)
                U(J,NM)=(Y*C)+(Z*S)
                U(J,I)=-(Y*S)+(Z*C)
  420         CONTINUE
            ENDIF
  430     CONTINUE
 2000     CONTINUE
          Z=W(K)
          IF(L.EQ.K)THEN
            IF(Z.LT.0.0D0)THEN
              W(K)=-Z
              DO 440 J=1,N
                V(J,K)=-V(J,K)
  440         CONTINUE
            ENDIF
            GOTO 3000
          ENDIF
          IF(ITS.EQ.30)THEN
            WRITE(6,'(1H ,A)')'NO CONVERGENCE IN 30 ITS.'
            RETURN
          ENDIF
          X=W(L)
          NM=K-1
          Y=W(NM)
          G=RV1(NM)
          H=RV1(K)
          F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0D0*H*Y)
          G=SQRT(F*F+1.0D0)
          F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X

          C=1.0D0
          S=1.0D0
          DO 470 J=L,NM
            I=J+1
            G=RV1(I)
            Y=W(I)
            H=S*G
            G=C*G
            Z=SQRT(F*F+H*H)
            RV1(J)=Z
            C=F/Z
            S=H/Z
            F= (X*C)+(G*S)
            G=-(X*S)+(G*C)
            H=Y*S
            Y=Y*C
            DO 450 JJ=1,N
              X=V(JJ,J)
              Z=V(JJ,I)
              V(JJ,J)= (X*C)+(Z*S)
              V(JJ,I)=-(X*S)+(Z*C)
  450       CONTINUE
            Z=SQRT(F*F+H*H)
            W(J)=Z
            IF(Z.NE.0.0D0)THEN
              Z=1.0D0/Z
              C=F*Z
              S=H*Z
            ENDIF
            F= (C*G)+(S*Y)
            X=-(S*G)+(C*Y)
            DO 460 JJ=1,M
              Y=U(JJ,J)
              Z=U(JJ,I)
              U(JJ,J)= (Y*C)+(Z*S)
              U(JJ,I)=-(Y*S)+(Z*C)
  460       CONTINUE
  470     CONTINUE
          RV1(L)=0.0D0
          RV1(K)=F
          W(K)=X
  480   CONTINUE
 3000   CONTINUE
  490 CONTINUE

      RETURN
      END
      SUBROUTINE SVBKSB(U,W,V,M,N,B,X)
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER MP,NP
      PARAMETER(MP=20,NP=MP)

      INTEGER I,J,K,M,N
      REAL*8   U(MP,NP),W(NP),V(NP,NP),B(MP),X(NP),TMP(100),S

      DO 100 J=1,N
        S=0.0D0
        IF(W(J).NE.0.0D0)THEN
          DO 110 I=1,M
            S=S+U(I,J)*B(I)
  110     CONTINUE
          S=S/W(J)
        ENDIF
        TMP(J)=S
  100 CONTINUE
      DO 200 J=1,N
        S=0.0D0
        DO 210 K=1,N
          S=S+V(J,K)*TMP(K)
  210   CONTINUE
        X(J)=S
  200 CONTINUE

      RETURN
      END
      SUBROUTINE JACOBI(AA,N,D,V,NROT,FLAG)
      IMPLICIT REAL*8(A-H,O-Z)
      INTEGER NP,NMAX
      PARAMETER(NP=20,NMAX=100)
      DIMENSION AA(NP,NP),A(NP,NP),D(NP),V(NP,NP),B(NMAX),Z(NMAX)
      CHARACTER*1 FLAG

      FLAG='0'
      DO 10 IP=1,N
        DO 20 IQ=1,N
          V(IP,IQ)=0.0D0
          A(IP,IQ)=AA(IP,IQ)
   20   CONTINUE
        V(IP,IP)=1.0D0
   10 CONTINUE
      DO 30 IP=1,N
        B(IP)=A(IP,IP)
        D(IP)=B(IP)
        Z(IP)=0.0D0
   30 CONTINUE
      NROT=0
      DO 100 I=1,50
        SM=0.0D0
        DO 110 IP=1,N-1
          DO 120 IQ=IP+1,N
            SM=SM+ABS(A(IP,IQ))
  120     CONTINUE
  110   CONTINUE
        IF(SM.EQ.0.0D0) RETURN
        IF(I.LT.4)THEN
          TRESH=0.2D0*SM/DFLOAT(N*N)
        ELSE
          TRESH=0.0D0
        ENDIF
        DO 200 IP=1,N-1
          DO 210 IQ=IP+1,N
            G=100.0D0*ABS(A(IP,IQ))
            IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP)))
     -                 .AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN
              A(IP,IQ)=0.0D0
            ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN
              H=D(IQ)-D(IP)
              IF(ABS(H)+G.EQ.ABS(H))THEN
                T=A(IP,IQ)/H
              ELSE
                THETA=0.5D0*H/A(IP,IQ)
                T=1.0D0/(ABS(THETA)+SQRT(1.0D0+THETA**2))
                IF(THETA.LT.0.0D0)T=-T
              ENDIF
              C=1.0D0/SQRT(1.0D0+T**2)
              S=T*C
              TAU=S/(1.0D0+C)
              H=T*A(IP,IQ)
              Z(IP)=Z(IP)-H
              Z(IQ)=Z(IQ)+H
              D(IP)=D(IP)-H
              D(IQ)=D(IQ)+H
              A(IP,IQ)=0.0D0
              DO 300 J=1,IP-1
                G=A(J,IP)
                H=A(J,IQ)
                A(J,IP)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
  300         CONTINUE
              DO 310 J=IP+1,IQ-1
                G=A(IP,J)
                H=A(J,IQ)
                A(IP,J)=G-S*(H+G*TAU)
                A(J,IQ)=H+S*(G-H*TAU)
  310         CONTINUE
              DO 320 J=IQ+1,N
                G=A(IP,J)
                H=A(IQ,J)
                A(IP,J)=G-S*(H+G*TAU)
                A(IQ,J)=H+S*(G-H*TAU)
  320         CONTINUE
              DO 330 J=1,N
                G=V(J,IP)
                H=V(J,IQ)
                V(J,IP)=G-S*(H+G*TAU)
                V(J,IQ)=H+S*(G-H*TAU)
  330         CONTINUE
              NROT=NROT+1
            ENDIF
  210     CONTINUE
  200   CONTINUE
        DO 400 IP=1,N
          B(IP)=B(IP)+Z(IP)
          D(IP)=B(IP)
          Z(IP)=0.0D0
  400   CONTINUE
  100 CONTINUE
      WRITE(6,*)'*** 50 iterations should never happen.'
      FLAG='1'
      RETURN
      END
      SUBROUTINE NABLAP(DPDPRM,PARM,W,H,V,XI)
      IMPLICIT REAL*8(A-Z)
*  FUNCTION
      EXTERNAL DQDU,DERU,ULIMIT
      REAL*8 DQDU

      INTEGER M
      PARAMETER(M=20)

      REAL*8 DPDPRM(1:4,M),PARM(M)
      REAL*8 W,H,V,XI
      INTEGER I,J,K,L
      REAL*8 GW(6),SIGMA
      REAL*8 U0,U1,U4
      REAL*8 DU0DG(M),DU1DG(M),DU4DG(M),DP0,DP1,DP4
      INTEGER ICON

      DO 120 I=1,4
        DO 110 J=1,M
          DPDPRM(I,J)=0.0D0
  110   CONTINUE
  120 CONTINUE
      DO 130 K=1,M
        DU0DG(K)=0.0D0
        DU1DG(K)=0.0D0
        DU4DG(K)=0.0D0
  130 CONTINUE

      GW(1)=-1.0D0
      GW(2)=PARM(1)
      GW(3)=PARM(2)
      GW(4)=PARM(3)
      GW(5)=PARM(4)
      GW(6)=PARM(5)
      SIGMA=PARM(6)

      CALL DERU(GW,SIGMA,W,H,V,XI,DU0DG,DU1DG,DU4DG)

      CALL ULIMIT(GW,SIGMA,W,H,V,XI,ICON,U0,U1,U4)

*	無業確率
*     MU(1)=Q(U0)
*	内職確率
*     MU(2)=Q(U1)-Q(U0)
*	雇用確率
*     MU(3)=Q(U4)-Q(U1)
*	雇用内職確率
*     MU(4)=1.0D0-Q(U4)
      DO 210 L=1,6
        DP0=DQDU(U0)*DU0DG(L+1)
        DP1=DQDU(U1)*DU1DG(L+1)
        DP4=DQDU(U4)*DU4DG(L+1)
        DPDPRM(1,L)= DP0
        DPDPRM(2,L)= DP1-DP0
        DPDPRM(3,L)= DP4-DP1
        DPDPRM(4,L)=    -DP4
  210 CONTINUE

      RETURN
      END
      SUBROUTINE DERU(GAMMA,SIGMA,W,H,V,XI,DFU0,DFU1,DFU4)
      IMPLICIT REAL*8(A-Z)
      INTEGER M
      PARAMETER(M=20)

      REAL*8 DFU0(M),DFU1(M),DFU4(M)
      REAL*8 GAMMA(6)
      REAL*8 G1,G2,G3,G4,G5,G4O,SIGMA
      REAL*8 W,H,V,XI
      REAL*8 XT
      SAVE
      DATA XT/1.0D0/

      E=DEXP(1.0D0)

      G1=GAMMA(1)
      G2=GAMMA(2)
      G3=GAMMA(3)
      G4=GAMMA(4)
      G5=GAMMA(5)
      G4O=GAMMA(6)

*     A:=G1*V*V-2*G3*V+G5;
*     B:=-(G1*V-G3)*XI-(G2+G3*XT)*V+G4O+G5*XT;
*     H0:=B/A;
*     H2:=G4/A*EXP((-1/2)*SIGMA**2);

*     U0:=(LOG(H0)-LOG(-H2))/SIGMA;

*     AA:=(G1*V*V-2*G3*V+G5)/2;
*     DD:=-(W-V)*(G1*V*V-2*G3*V+G5)*H*(G1*((W+V)*H+2*XI)
*         +2*G2+2*G3*(XT-H));
*     HDQ1:=H+SQRT(DD)/(2*AA);
*     U1:=(LOG(H0-HDQ1)-LOG(-H2))/SIGMA;

*     AAA:=G1*V*V-2*G3*V+G5;
*     BBB:=(G1*V-G3)*(W-V)*H;
*     HDQ4:=H+BBB/AAA;
*     U4:=(LOG(H0-HDQ4)-LOG(-H2))/SIGMA;

*     DF(U0,G2);
*     DF(U0,G3);
*     DF(U0,G4);
*     DF(U0,G5);
*     DF(U0,G4O);
*     DF(U0,SIGMA);

*     DF(U1,G2);
*     DF(U1,G3);
*     DF(U1,G4);
*     DF(U1,G5);
*     DF(U1,G4O);
*     DF(U1,SIGMA);

*     DF(U4,G2);
*     DF(U4,G3);
*     DF(U4,G4);
*     DF(U4,G5);
*     DF(U4,G4O);
*     DF(U4,SIGMA);


*      DF(U0,G2);
      DFU0(2)=V/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5-G4O+V*G2))

*      DF(U0,G3);
      DFU0(3)=-(XI-XT*V)/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5-
     . G4O+V*G2))

*      DF(U0,G4);
      DFU0(4)=-1./(SIGMA*G4)

*      DF(U0,G5);
      DFU0(5)=-XT/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5-G4O+V*G2
     . ))

*      DF(U0,G4O);
      DFU0(6)=-1./(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5-G4O+V*G2
     . ))

*      DF(U0,SIGMA);
      DFU0(7)=(LOG((-G4)/(SQRT(E)**(SIGMA**2)*V**2*G1-2.*SQRT(
     . E)**(SIGMA**2)*V*G3+SQRT(E)**(SIGMA**2)*G5))-LOG((-
     . XI*V*G1+XI*G3-XT*V*G3+XT*G5+G4O-V*G2)/(V**2*G1-2.*V
     . *G3+G5))+SIGMA**2)/SIGMA**2


*      DF(U1,G2);
      ANS1=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*V+H*V**3*G1-H*V**2*W*G1-2.*H*V**2*G3+2.*H
     . *V*W*G3+H*V*G5-H*W*G5
      ANS12=2.*H**2*V**2*W*G1*G3+H**2*V**2*G1*G5+4.*H**2*V
     . **2*G3**2+2.*H**2*V*W**2*G1*G3-4.*H**2*V*W*G3**2-2.
     . *H**2*V*G3*G5-H**2*W**2*G1*G5+2.*H**2*W*G3*G5+2.*H*
     . V**3*G1*G2-2.*H*V**2*W*G1*G2-4.*H*V**2*G3*G2+4.*H*V
     . *W*G3*G2+2.*H*V*G2*G5-2.*H*W*G2*G5
      ANS11=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*V*G2+2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1
     . **2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*
     . G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V
     . **2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.
     . *XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H
     . **2*V**3*G1*G3-H**2*V**2*W**2*G1**2+ANS12
      ANS10=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*G4O+ANS11
      ANS9=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*G5+ANS10
      ANS8=-2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**
     . 2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1
     . *G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2
     . *W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT
     . *H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*
     . V**3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*
     . G3+H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**
     . 2*G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2
     . *G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*
     . G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-
     . 2.*H*W*G2*G5)*H*V*G3+ANS9
      ANS7=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*V**2*G1+ANS8
      ANS6=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*G5+ANS7
      ANS5=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*V*G3+ANS6
      ANS4=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*G3+ANS5
      ANS3=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*V*G1+ANS4
      ANS2=SIGMA*ANS3
      DFU1(2)=ANS1/ANS2

*      DF(U1,G3);
      ANS4=2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2
     . -4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*
     . G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*
     . W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*
     . H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V
     . **3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3
     . +H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*
     . G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*
     . G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*
     . G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-
     . 2.*H*W*G2*G5)*H*V+2.*XI*H*V**2*G1-2.*XI*H*V*W*G1-XT*
     . H*V**3*G1+XT*H*V**2*W*G1+4.*XT*H*V**2*G3-4.*XT*H*V*
     . W*G3-XT*H*V*G5+XT*H*W*G5+2.*H**2*V**3*G1-H**2*V**2*
     . W*G1-4.*H**2*V**2*G3-H**2*V*W**2*G1+4.*H**2*V*W*G3+
     . H**2*V*G5-H**2*W*G5+2.*H*V**2*G2-2.*H*V*W*G2
      ANS3=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*V+ANS4
      ANS2=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI+ANS3
      ANS15=2.*H**2*V**2*W*G1*G3+H**2*V**2*G1*G5+4.*H**2*V
     . **2*G3**2+2.*H**2*V*W**2*G1*G3-4.*H**2*V*W*G3**2-2.
     . *H**2*V*G3*G5-H**2*W**2*G1*G5+2.*H**2*W*G3*G5+2.*H*
     . V**3*G1*G2-2.*H*V**2*W*G1*G2-4.*H*V**2*G3*G2+4.*H*V
     . *W*G3*G2+2.*H*V*G2*G5-2.*H*W*G2*G5
      ANS14=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*V*G2+2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1
     . **2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*
     . G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V
     . **2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.
     . *XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H
     . **2*V**3*G1*G3-H**2*V**2*W**2*G1**2+ANS15
      ANS13=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*G4O+ANS14
      ANS12=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*G5+ANS13
      ANS11=-2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1
     . **2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*
     . G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V
     . **2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.
     . *XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H
     . **2*V**3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*
     . G1*G3+H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*
     . W**2*G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W
     . **2*G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2
     . *W*G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*
     . G5-2.*H*W*G2*G5)*H*V*G3+ANS12
      ANS10=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*V**2*G1+ANS11
      ANS9=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*G5+ANS10
      ANS8=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*V*G3+ANS9
      ANS7=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*G3+ANS8
      ANS6=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*V*G1+ANS7
      ANS5=SIGMA*ANS6
      ANS1=ANS2/ANS5
      DFU1(3)=-ANS1

*      DF(U1,G4);
      DFU1(4)=-1./(SIGMA*G4)

*      DF(U1,G5);
      ANS3=-2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**
     . 2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1
     . *G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2
     . *W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT
     . *H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*
     . V**3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*
     . G3+H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**
     . 2*G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2
     . *G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*
     . G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-
     . 2.*H*W*G2*G5)*H-2.*XI*H*V*G1+2.*XI*H*W*G1-2.*XT*H*V*
     . G3+2.*XT*H*W*G3-H**2*V**2*G1+2.*H**2*V*G3+H**2*W**2
     . *G1-2.*H**2*W*G3-2.*H*V*G2+2.*H*W*G2
      ANS2=2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2
     . -4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*
     . G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*
     . W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*
     . H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V
     . **3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3
     . +H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*
     . G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*
     . G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*
     . G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-
     . 2.*H*W*G2*G5)*XT+ANS3
      ANS14=2.*H**2*V**2*W*G1*G3+H**2*V**2*G1*G5+4.*H**2*V
     . **2*G3**2+2.*H**2*V*W**2*G1*G3-4.*H**2*V*W*G3**2-2.
     . *H**2*V*G3*G5-H**2*W**2*G1*G5+2.*H**2*W*G3*G5+2.*H*
     . V**3*G1*G2-2.*H*V**2*W*G1*G2-4.*H*V**2*G3*G2+4.*H*V
     . *W*G3*G2+2.*H*V*G2*G5-2.*H*W*G2*G5
      ANS13=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*V*G2+2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1
     . **2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*
     . G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V
     . **2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.
     . *XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H
     . **2*V**3*G1*G3-H**2*V**2*W**2*G1**2+ANS14
      ANS12=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*G4O+ANS13
      ANS11=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*G5+ANS12
      ANS10=-2.*SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1
     . **2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*
     . G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V
     . **2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.
     . *XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H
     . **2*V**3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*
     . G1*G3+H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*
     . W**2*G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W
     . **2*G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2
     . *W*G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*
     . G5-2.*H*W*G2*G5)*H*V*G3+ANS11
      ANS9=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*H*V**2*G1+ANS10
      ANS8=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*G5+ANS9
      ANS7=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XT*V*G3+ANS8
      ANS6=-SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-
     . 4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5
     . -2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*
     . G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*
     . V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**
     . 3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H
     . **2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1
     . *G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*G3+ANS7
      ANS5=SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.
     . *XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-
     . 2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1
     . *G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)*XI*V*G1+ANS6
      ANS4=2.*SIGMA*ANS5
      ANS1=ANS2/ANS4
      DFU1(5)=-ANS1

*      DF(U1,G4O);
      DFU1(6)=-1./(SIGMA*(SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2
     . *W*G1**2-4.*XI*H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI
     . *H*V*G1*G5-2.*XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT
     . *H*V**2*W*G1*G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**
     . 2+2.*XT*H*V*G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-
     . 4.*H**2*V**3*G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2
     . *W*G1*G3+H**2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2
     . *V*W**2*G1*G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**
     . 2*W**2*G1*G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V
     . **2*W*G1*G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*
     . G2*G5-2.*H*W*G2*G5)+XI*V*G1-XI*G3+XT*V*G3-XT*G5+H*V
     . **2*G1-2.*H*V*G3+H*G5-G4O+V*G2))

*      DF(U1,SIGMA);
      DFU1(7)=(LOG((-G4)/(SQRT(E)**(SIGMA**2)*V**2*G1-2.*SQRT(
     . E)**(SIGMA**2)*V*G3+SQRT(E)**(SIGMA**2)*G5))-LOG((-
     . SQRT(2.*XI*H*V**3*G1**2-2.*XI*H*V**2*W*G1**2-4.*XI*
     . H*V**2*G1*G3+4.*XI*H*V*W*G1*G3+2.*XI*H*V*G1*G5-2.*
     . XI*H*W*G1*G5+2.*XT*H*V**3*G1*G3-2.*XT*H*V**2*W*G1*
     . G3-4.*XT*H*V**2*G3**2+4.*XT*H*V*W*G3**2+2.*XT*H*V*
     . G3*G5-2.*XT*H*W*G3*G5+H**2*V**4*G1**2-4.*H**2*V**3*
     . G1*G3-H**2*V**2*W**2*G1**2+2.*H**2*V**2*W*G1*G3+H**
     . 2*V**2*G1*G5+4.*H**2*V**2*G3**2+2.*H**2*V*W**2*G1*
     . G3-4.*H**2*V*W*G3**2-2.*H**2*V*G3*G5-H**2*W**2*G1*
     . G5+2.*H**2*W*G3*G5+2.*H*V**3*G1*G2-2.*H*V**2*W*G1*
     . G2-4.*H*V**2*G3*G2+4.*H*V*W*G3*G2+2.*H*V*G2*G5-2.*H
     . *W*G2*G5)-XI*V*G1+XI*G3-XT*V*G3+XT*G5-H*V**2*G1+2.*
     . H*V*G3-H*G5+G4O-V*G2)/(V**2*G1-2.*V*G3+G5))+SIGMA**
     . 2)/SIGMA**2


*      DF(U4,G2);
      DFU4(2)=V/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5+H*V*W*G1-H
     . *V*G3-H*W*G3+H*G5-G4O+V*G2))

*      DF(U4,G3);
      DFU4(3)=-(XI-XT*V+H*V+H*W)/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3
     . -XT*G5+H*V*W*G1-H*V*G3-H*W*G3+H*G5-G4O+V*G2))

*      DF(U4,G4);
      DFU4(4)=-1./(SIGMA*G4)

*      DF(U4,G5);
      DFU4(5)=-(XT-H)/(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5+H*V*
     . W*G1-H*V*G3-H*W*G3+H*G5-G4O+V*G2))

*      DF(U4,G4O);
      DFU4(6)=-1./(SIGMA*(XI*V*G1-XI*G3+XT*V*G3-XT*G5+H*V*W*G1
     . -H*V*G3-H*W*G3+H*G5-G4O+V*G2))

*      DF(U4,SIGMA);
      DFU4(7)=(LOG((-G4)/(SQRT(E)**(SIGMA**2)*V**2*G1-2.*SQRT(
     . E)**(SIGMA**2)*V*G3+SQRT(E)**(SIGMA**2)*G5))-LOG((-
     . XI*V*G1+XI*G3-XT*V*G3+XT*G5-H*V*W*G1+H*V*G3+H*W*G3-
     . H*G5+G4O-V*G2)/(V**2*G1-2.*V*G3+G5))+SIGMA**2)/
     . SIGMA**2


      RETURN
      END
      FUNCTION DQDU(U)
      IMPLICIT REAL*8(A-Z)

      REAL*8 U,PI,DQDU
      COMMON/BLKPI/PI
      SAVE

      IF(ABS(U).GE.10.0)THEN
        DQDU=0.0
      ELSE
        DQDU=-1.0D0/SQRT(2.0D0*PI)*DEXP(-U*U/2.0D0)
      ENDIF

      RETURN
      END
      SUBROUTINE LS4(MU,GAMMA,SIGMA,W,H,V,XI,ICON,U0,U1,U4)
      IMPLICIT REAL*8(A-Z)

      EXTERNAL ULIMIT

      REAL*8 MU(4)
      REAL*8 GAMMA(6),SIGMA
      REAL*8 W,H,V,XI
      INTEGER ICON
      REAL*8 U0,U1,U4

      CALL ULIMIT(GAMMA,SIGMA,W,H,V,XI,ICON,U0,U1,U4)
      IF(ICON.NE.0)RETURN

*	無業確率
      MU(1)=Q(U0)
*	内職確率
      MU(2)=Q(U1)-MU(1)
*	雇用確率
      MU(3)=Q(U4)-MU(2)-MU(1)
*	雇用内職確率
      MU(4)=1.0D0-Q(U4)

      RETURN
      END
      SUBROUTINE ULIMIT(GAMMA,SIGMA,W,H,V,XI,ICON,U0,U1,U4)
      IMPLICIT REAL*8(A-Z)

      EXTERNAL HDQ1,HDQ4
      REAL*8 HDQ1,HDQ4

      INTEGER ICON
      REAL*8 W,H,V,XI
      REAL*8 H0,H2
      REAL*8 G1,G2,G3,G4BAR,G4O,SIGMA
      REAL*8 U0,U1,U4
      REAL*8 GAMMA(6)
      REAL*8 HDQ1X,HDQ4X
      REAL*8 T
      REAL*8 A,B
      SAVE
      DATA T/1.0D0/

      G1   =GAMMA(1)
      G2   =GAMMA(2)
      G3   =GAMMA(3)
      G4BAR=GAMMA(4)
      G5   =GAMMA(5)
      G4O  =GAMMA(6)

      A=G1*V*V-2.0D0*G3*V+G5
      B=-(G1*V-G3)*XI-(G2+G3*T)*V+G4O+G5*T
      H0=B/A
      H2=G4BAR/A*EXP(-0.5D0*SIGMA**2)

      U0=0.0D0
      U1=0.0D0
      U4=0.0D0

      U0=(LOG(H0)-LOG(-H2))/SIGMA

      HDQ1X=HDQ1(GAMMA,W,H,V,XI,ICON)
      IF(ICON.NE.0)RETURN
      U1=(LOG(H0-HDQ1X)-LOG(-H2))/SIGMA

      HDQ4X=HDQ4(GAMMA,W,H,V)
      U4=(LOG(H0-HDQ4X)-LOG(-H2))/SIGMA

      RETURN
      END
      FUNCTION HDQ1(GAMMA,W,H,V,XI,ICON)
      IMPLICIT REAL*8(A-Z)
      REAL*8 HDQ1
      REAL*8 GAMMA(6)
      REAL*8 W,H,V,XI
      INTEGER ICON
      REAL*8 G1,G2,G3,G4BAR,G5,G4O
      REAL*8 T
      REAL*8 A,D
      SAVE
      DATA T/1.0D0/

      G1   =GAMMA(1)
      G2   =GAMMA(2)
      G3   =GAMMA(3)
      G4BAR=GAMMA(4)
      G5   =GAMMA(5)
      G4O  =GAMMA(6)

      A=(G1*V*V-2.0D0*G3*V+G5)/2.0D0
*      B=-(G1*V*V-2.0D0*G3*V+G5)*H
*      C=(G1*W*W-2.0D0*G3*W+G5)*H*H/2.0D0+(W-V)*(G1*XI+G2+G3*T)*H
*      D=B*B-4*A*C
      D=-(W-V)*(G1*V*V-2.0D0*G3*V+G5)*H*(G1*((W+V)*H+2.0D0*XI)
     :   +2.0D0*G2+2.0D0*G3*(T-H))
      ICON=0
      IF(D.LT.0.0D0)THEN
        ICON=1
        HDQ1=0.0
      ELSE
        HDQ1=H+SQRT(D)/(2.0D0*A)
        IF(HDQ1.GT.0.0D0.AND.HDQ1.LT.H)THEN
        ELSE
          HDQ1=0.0D0
          ICON=2
        ENDIF
      ENDIF

      RETURN
      END
      FUNCTION HDQ4(GAMMA,W,H,V)
      IMPLICIT REAL*8(A-Z)
      REAL*8 HDQ4
      REAL*8 GAMMA(6)
      REAL*8 W,H,V
      REAL*8 G1,G2,G3,G4BAR,G5,G4O
      REAL*8 T
      REAL*8 A,B
      SAVE
      DATA T/1.0D0/

      G1   =GAMMA(1)
      G2   =GAMMA(2)
      G3   =GAMMA(3)
      G4BAR=GAMMA(4)
      G5   =GAMMA(5)
      G4O  =GAMMA(6)

      A=G1*V*V-2.0D0*G3*V+G5
      B=(G1*V-G3)*(W-V)*H
      HDQ4=H+B/A

      RETURN
      END
      FUNCTION Q(U)
      IMPLICIT REAL*8(A-Z)

*     /*  標準正規分布の上側確率     */
*     /*   P(U>U*)                   */
*     /*                             */

      REAL*8   Q,U,V,V2
      REAL*8   P,SUM,THETA(0:3),TERM,TERMP
      REAL*8   PI,EPS
      INTEGER  R
      CHARACTER*1 FLAG

      COMMON/BLKEPS/EPS
      COMMON/BLKPI/PI
      SAVE

      IF(ABS(U).GT.10.0D0)THEN
        IF(U.GT.0.0D0)THEN
          Q=0.0D0
        ELSE
          Q=1.0D0
        ENDIF
        RETURN
      ENDIF

      V =U/2.0D0
      V2=V*V
      THETA(0)=1.0D0
      THETA(1)=V2
      SUM=1.0D0
      TERMP=1.0D0
      R=0
  100 CONTINUE
        R=R+1
        THETA(2)=V2*(THETA(1)-THETA(0))/DFLOAT(2*R  )
        THETA(3)=V2*(THETA(2)-THETA(1))/DFLOAT(2*R+1)
        TERM=THETA(2)/DFLOAT(2*R+1)
        SUM=SUM+TERM
        FLAG='0'
        IF(ABS(TERMP).LT.EPS)FLAG='1'
        TERMP=TERM
        THETA(0)=THETA(2)
        THETA(1)=THETA(3)
      IF(.NOT.(ABS(TERM).LT.EPS.AND.FLAG.EQ.'1'))GOTO 100

      P=U/SQRT(2.0D0*PI)*EXP((-1.0D0)*U**2/8.0D0)*SUM

      Q=0.5D0-P

      RETURN
      END
      INTEGER FUNCTION PRECHK(PARM,LVRP,IDTR,KT)
      IMPLICIT REAL*8(A-Z)
      INTEGER L,M
      PARAMETER(L=200,M=20)

* FUNCTIONS
      INTEGER THCHK4

      INTEGER KR
      PARAMETER(KR=70)
      INTEGER KT,IDTR(KR)
      REAL*8 CRTR(KR)

      REAL*8 PARM(M),GW(6),SGMW
      INTEGER LVRP,CODEW,CODE,I
      INTEGER NPRM,NPRMS,LPRMS(M)
      INTEGER NSMPL
      REAL*8 WW,HW,VW,XI
      REAL*8 WWTAB(L)
      REAL*8 HWTAB(L)
      REAL*8 VWTAB(L)
      REAL*8 XITAB(L)
      REAL*8 WEITAB(L)
      COMMON/BLKDT/WWTAB,HWTAB,VWTAB,XITAB,WEITAB
      COMMON/BLKCL/NPRM,NPRMS,LPRMS
      COMMON/BLKMT/NSMPL
      SAVE

      GW(1)=-1.0D0
      GW(2)=PARM(1)
      GW(3)=PARM(2)
      GW(4)=PARM(3)
      GW(5)=PARM(4)
      GW(6)=PARM(5)
      SGMW =PARM(6)

      CODE=0
      DO 100 I=1,NSMPL
        IF(CODE.EQ.1)GOTO 110
        IF(LVRP.GT.1)WRITE(6,'(1H ,A,I5)')'========== GROUP',I
        WW=WWTAB(I)
        HW=HWTAB(I)
        VW=VWTAB(I)
        XI=XITAB(I)
        CODEW=THCHK4(LVRP,GW,SGMW,WW,HW,VW,XI,IDTR,CRTR,KT)
        IF(CODEW.EQ.1)CODE=1
  100 CONTINUE
  110 CONTINUE

      PRECHK=CODE

      RETURN
      END
      INTEGER FUNCTION THCHK4(LVRP,GAMMA,SIGMA,W,H,V,XI,IDTR,CRTR,KT)

********************************************************
***
***		4者択一モデルの理論制約
***
***		理論制約を充足する場合には THCHK4=0 の値を返す
***
********************************************************

      IMPLICIT REAL*8(A-Z)
      INTEGER LVRP
      REAL*8 GAMMA(6),SIGMA
      REAL*8 W,H,V,XI
      REAL*8 T
      REAL*8 G1,G2,G3,G4,G5,G4O
      REAL*8 WDUM,HDUM,VDUM,XIDUM
      REAL*8 EPSMAX,EPSMIN,LMDMAX,LMDMIN
      REAL*8 DTRMNT,EEPS,UM,RM,DUDLMD,DUDX
      REAL*8 CENTR1,CENTR2,ALPHA
      REAL*8 H0,H2,ULMT
      CHARACTER*10 SHAPE(3)
      CHARACTER*10 CURVE
      CHARACTER*1  FAIL0,FAIL1,FAIL2,FAIL3,FAIL4,FAILU

      INTEGER KR
      PARAMETER(KR=70)
      INTEGER KT,IDTR(KR)
      REAL*8 CRTR(KR)

      REAL*8 FDUDL,FDUDX
      REAL*8 FCNVEX
      REAL*8 EPSLON,LAMBDA
      REAL*8 DUDL,DUDXI


      INTEGER  I,ISGM,ITHRS,MU,MR,CODE
      SAVE
      DATA T/1.0D0/


      FDUDL(EPSLON,LAMBDA) = G3*XIDUM+G4O+G4*EPSLON+G5*LAMBDA
      FDUDX(LAMBDA)        = G1*XIDUM+G2+G3*LAMBDA
      FCNVEX(DUDL,DUDXI)    = 2.0D0*G3*DUDL*DUDXI-G5*DUDXI**2-G1*DUDL**2


      SHAPE(1)='ELLIPSE'
      SHAPE(2)='HYPERBOLA'
      SHAPE(3)='LINE X'

      LMDMIN=0.25D0*T
      LMDMAX=1.0D0 *T

      G1=GAMMA(1)
      G2=GAMMA(2)
      G3=GAMMA(3)
      G4=GAMMA(4)
      G5=GAMMA(5)
      G4O=GAMMA(6)

      KT=0

*  /*     CHECK THE THEORETICAL RESTRICTIONS 1=A1<0,
*         1-2=SLOPE RIGHT DOWN,3=MU OF INCOME,4-6=MU OF LEISURE,
*         7-18=CONVEX TO ORIGINE,19=H0>IMAX,20=H2<0    */
      FAIL0='0'
      FAILU='0'

      IF(SIGMA.LE.0.0D0.OR.SIGMA.GT.10.0D0)THEN
        KT=KT+1
        IDTR(KT)=99
        IF(LVRP.GT.1)WRITE(6,'(1H ,A,G12.6)')'XXXX WARING:SIGMA=',SIGMA
        THCHK4=1
        RETURN
      ENDIF

*  /* REQUIRED CONDITION FOR BEING ELLIPSE   : G1*G5-G3*G3>0
*                                  HYPERBOLA : G1*G5-G3*G3<0 */

*	無差別曲線を双曲線に限定する
      DTRMNT=G1*G5-G3**2
      IF(DTRMNT.GE.0.0D0)THEN
        FAILU='1'
        THCHK4=1
        RETURN
      ENDIF
      IF(DTRMNT.GT.0.0D0)THEN
        CURVE=SHAPE(1)
      ELSE IF(DTRMNT.LT.0.0D0)THEN
        CURVE=SHAPE(2)
      ELSE
        CURVE=SHAPE(3)
      ENDIF

      IF(LVRP.GT.1)THEN
        WRITE(6,'(1H ,4(A,I1,A,G13.6)/1H ,2(A,G13.6))')
     :           (' G',I,'=',GAMMA(I),I=2,5),
     :            ' G4O=',GAMMA(6),' SIGMA=',SIGMA
        WRITE(6,'(1H ,13X,A)')'|G1 G3|'
        WRITE(6,'(1H ,A,G15.7,A,A)')
     :          ' DETERMINANT |G3 G5| =',DTRMNT,
     :          '    INDIFERENCE CURVE = ',CURVE
        IF(FAIL0.EQ.'1')
     :    WRITE(6,'(1H ,A)')'XXXX WARNING:INDIF. CURVE NOT ELLIPSE.'
        CENTR1=(G3*(G4O+G4*EXP(-0.5D0*SIGMA**2))-G5*G2)/DTRMNT
        CENTR2=(G3*G2-G1*(G4O+G4*EXP(-0.5D0*SIGMA**2)))/DTRMNT
        ALPHA=ATAN(2*G3/(G1-G5))/2.0/3.141593*180.0
        WRITE(6,'(1H ,A,G15.8,A,G15.8,A/1H ,5X,A,F10.5,A)')
     :                'ﾁﾕｳｼﾝ ﾉ ｻﾞﾋﾖｳ = (',CENTR1,',',CENTR2,
     :     ') WHEN U*=0.0,','ｼﾞｸ ﾉ X-ｼﾞｸ ﾆ ﾀｲｽﾙ ｶｸﾄﾞ = ',ALPHA,'ﾄﾞ'
      ENDIF

      XIDUM=XI
      HDUM=H
      WDUM=W
      VDUM=V
        FAIL1='0'
        FAIL2='0'
        FAIL3='0'
        IF(LVRP.GT.1)THEN
          WRITE(6,'(1H ,A,F10.6,2(A,F8.2))')
     :         'HOUR ASSIGNED=',HDUM,' WAGE RATE=',WDUM,
     :         ' PRINCIPAL INCOME=',XIDUM
          WRITE(6,'(1H ,A)')'THEORETICAL RESTRICTIONS : 1=A1<0,',
     :     '1-2=SLOPE RIGHT DOWN,3=MU OF INCOME,4-6=MU OF LEISURE,',
     :     '7-18=CONVEX TO ORIGINE,19=H0>IMAX,20=H2<0'
          WRITE(6,'(1H ,A)')'4者択一モデルの理論制約',
     :     '21=φ関数の勾配正,22-23=0<φ[H(d)=0]=f[H(d)=0]',
     :     '24=φ[H(d)=0]=f[H(d)=0]<h,25=q1の座標<h,26=0<q1の座標',
     :     '27=q4の座標<a,28=h<q4の座標,29=q1の計算の平方根,30-31=G4>0'
        ENDIF


***
***		2者択一モデルの理論制約
***

        DO 10 I=1,40
          CRTR(I)=1.0D0
   10   CONTINUE

        CRTR(1)=WDUM-G3/G1
        CRTR(2)=WDUM**2-2.0D0*WDUM*G3/G1+G5/G1

        IF(G3.GE.0.0D0)THEN
          CRTR(3)=G1*XIDUM+G2+G3*LMDMIN
        ELSE
          CRTR(3)=G1*XIDUM+G2+G3*LMDMAX
        ENDIF

        EEPS=EXP(-0.5D0*SIGMA**2)

        DO 200 ISGM=1,2
          EPSMAX=EEPS*EXP(DFLOAT( ISGM)*SIGMA)
          EPSMIN=EEPS*EXP(DFLOAT(-ISGM)*SIGMA)
          IF(G5.GT.0.0D0)THEN
            IF(G4.GT.0.0D0)THEN
              CRTR(3+ISGM)=G3*XIDUM+G4O+G4*EPSMIN+G5*LMDMIN
            ELSE
              CRTR(3+ISGM)=G3*XIDUM+G4O+G4*EPSMAX+G5*LMDMIN
            ENDIF
          ELSE
            IF(G4.GT.0.0D0)THEN
              CRTR(3+ISGM)=G3*XIDUM+G4O+G4*EPSMIN+G5*LMDMAX
            ELSE
              CRTR(3+ISGM)=G3*XIDUM+G4O+G4*EPSMAX+G5*LMDMAX
            ENDIF
          ENDIF
  200   CONTINUE

        ITHRS=6
        DO 300 ISGM=1,3

          EPSMAX=EEPS*EXP(DFLOAT( ISGM)*SIGMA)
          EPSMIN=EEPS*EXP(DFLOAT(-ISGM)*SIGMA)

          DUDLMD = FDUDL(EPSMAX,LMDMAX)
          DUDX   = FDUDX(LMDMAX)
          ITHRS=ITHRS+1
          CRTR(ITHRS) = FCNVEX(DUDLMD,DUDX)

          DUDLMD = FDUDL(EPSMAX,LMDMIN)
          DUDX   = FDUDX(LMDMIN)
          ITHRS=ITHRS+1
          CRTR(ITHRS) = FCNVEX(DUDLMD,DUDX)

          DUDLMD = FDUDL(EPSMIN,LMDMAX)
          DUDX   = FDUDX(LMDMAX)
          ITHRS=ITHRS+1
          CRTR(ITHRS) = FCNVEX(DUDLMD,DUDX)

          DUDLMD = FDUDL(EPSMIN,LMDMIN)
          DUDX   = FDUDX(LMDMIN)
          ITHRS=ITHRS+1
          CRTR(ITHRS) = FCNVEX(DUDLMD,DUDX)

  300   CONTINUE

        CRTR(19)=(G4O-G2*WDUM-G3*WDUM*(T-HDUM)+G5*(T-HDUM/2.0D0)
     :          -G1*WDUM**2*HDUM/2.0D0)/(G1*WDUM-G3)-XIDUM
        CRTR(20)=-G4/(G1*WDUM-G3)

        DO 600 ITHRS=1,20
          IF(CRTR(ITHRS).LE.0.0D0)THEN
            IF(ITHRS.LE.6)THEN
              FAIL1='1'
            ELSE IF(ITHRS.LE.18)THEN
              FAIL2='1'
            ELSE
              FAIL3='1'
            ENDIF
            KT=KT+1
            IDTR(KT)=ITHRS
          ENDIF
  600   CONTINUE

        IF(FAIL3.EQ.'0')THEN
          IF(LVRP.GT.1)THEN
            H0=CRTR(19)+XIDUM
            H2=-CRTR(20)*EXP(-0.5D0*SIGMA**2)
            ULMT=(LOG(H0-XIDUM)-LOG(-H2))/SIGMA
            WRITE(6,'(1H ,2(1X,A,G13.6),1X,A,G13.6)')
     :               'H0=',H0,'H2=',H2,'U=',ULMT
          ENDIF
        ENDIF

        IF(FAIL1.EQ.'1'.OR.FAIL2.EQ.'1'.OR.FAIL3.EQ.'1')FAILU='1'


***
***		4者択一モデルの理論制約
***

*	φ関数の勾配正
        CRTR(21)=(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
     :          *(G1*WDUM*WDUM-2.0D0*G3*WDUM+G5)
*	0<φ[H(d)=0]=f[H(d)=0]
        CRTR(22)=-(G1*WDUM*WDUM-2.0D0*G3*WDUM+G5)
        CRTR(23)=-(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
*	φ[H(d)=0]=f[H(d)=0]<h
        CRTR(24)=-0.5D0*HDUM*(G1*WDUM*WDUM-2.0D0*G3*WDUM+G5)
     :          -(WDUM-VDUM)*(G1*XIDUM+G2+G3*LMDMAX)
*	q1の座標<h
        CRTR(25)=-HDUM*HDUM*(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
     :          +HDUM*((G1*WDUM*WDUM-2.0D0*G3*WDUM+G5)*HDUM
     :          +2.0D0*(WDUM-VDUM)*(G1*XIDUM+G2+G3*LMDMAX))
*	0<q1の座標
        CRTR(26)=-HDUM*((G1*WDUM*WDUM-2.0D0*G3*WDUM+G5)*HDUM
     :          +2.0D0*(WDUM-VDUM)*(G1*XIDUM+G2+G3*LMDMAX))
*	q4の座標<a
*        HDMAX=(G5*LMDMAX-(G1*VDUM-G3)*XIDUM-VDUM*(G2+G3*LMDMAX))  (G4O=0のケース)
*     :       /(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
        HDMAX=(G4O+G5*LMDMAX-(G1*VDUM-G3)*XIDUM-VDUM*(G2+G3*LMDMAX))
     :       /(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
        HDQ4=HDUM+(G1*VDUM-G3)*(WDUM-VDUM)*HDUM
     :      /(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)
        CRTR(27)=HDMAX-HDQ4
*	h<q4の座標
        CRTR(28)=-(G1*VDUM-G3)
*	q1の計算の平方根
        CRTR(29)=-(WDUM-VDUM)*(G1*VDUM*VDUM-2.0D0*G3*VDUM+G5)*HDUM
     :   *(G1*((WDUM+VDUM)*HDUM+2.0D0*XIDUM)+2.0D0*G2+2.0D0*G3*(T-HDUM))
*	G4>0
        CRTR(30)=G4
*	G4O>0
        CRTR(31)=G4O


        FAIL4='0'
        DO 700 ITHRS=21,31
          IF(CRTR(ITHRS).LE.0.0D0)THEN
            FAIL4='1'
            KT=KT+1
            IDTR(KT)=ITHRS
          ENDIF
  700   CONTINUE

        IF(LVRP.GT.1)THEN
          WRITE(6,'(1H ,A,28I3)') 'VIOLATE THE RESTRICTION:',
     :             (IDTR(ITHRS),ITHRS=1,KT)
          WRITE(6,'(1H ,A,A)')'THEORETICAL RESTRICTIONS : R.I >0 ',
     :                        'FOR ALL I=1,2,3, ...,20 & 21,...,31'
          WRITE(6,'(4(1X,A2,I2,A1,G13.6))')
     :             ('R.',I,':',CRTR(I),I=1,31)
        ENDIF

        IF(FAIL4.EQ.'1')FAILU='1'

      IF(FAILU.EQ.'1')THEN
        CODE=1
      ELSE
        CODE=0
        IF(LVRP.GT.1) WRITE(6,'(1H ,A)')
     :  '****** PARAMETERS PASSED ALL THE THEORETICAL RESTRICTIONS.'
      ENDIF

      THCHK4=CODE

      RETURN
      END
