      program ssmlmnp
c     ported to the sun by vah, Feb 91
c     program statement also added
c
c     vah 18:55:22.11, Thu, Mar 7, 1991
c
c     1. lower case
c     2. fix '$include' statements
c     3. replace suitably 'call getarg' to 'call parsecl'
c     4. correct '\' in format statements
c     5. comment out the '$large' statement
c
c-----sml.for-----------------------------------------------------------
c
c     =====================================================
c     ****    multinomial multiperiod probit model     ****
c     * estimation by smooth simulated maximum likelihood *
c     =====================================================
c
c     (c) axel boersch-supan    version 2.1   1. september 1990
c
c     based on: - geweke, 1989
c               - hajivassiliou and mcfadden, 1990
c.pl100/.rm72/
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      real*4      data(maxbuf)
      integer*4   iseed
      character*8 pname(maxnpp)
      character*20 runtag
c
c     the following list includes all named common blocks
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / dname / pname
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / doptm / mchk0,mchk1,iter1,accu1,ialg1,
     *                       mchk2,iter2,accu2,ialg2,nwesml
      common / dsimu / nsimul
      common / dseed / iseed
      common / dtag  / runtag
      common / ddata / data
      common / dfreq / wtsum,share(maxnt),wshare(maxnt),
     *                       amean(maxndm),wmean(maxndm)
      common / dxvar / x(maxnt,maxnpp),ic(maxper),wt
      common / ddiff / idiff2
      common / dival / ival
      common / doutp / kanal0,ispool,ipause
c
c     reduce stack size
c
      common / dfunc / em1(maxnt1),em2(maxnt1),sg1(maxnt1,maxnt1),
     *                 dem1(maxnt1,maxnp),dem2(maxnt1,maxnp),
     *                 dsg1(maxnt1,maxnt1,maxnp),da(maxnp),dalf(maxnp)
      common / dzbar / zquer(maxnt1),covz(maxnt1,maxnt1),
     *                 dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns)
      common / dopt1 / parm(maxnp),grad(maxnp),hess(maxnp,maxnp),
     *                 hinv(maxnp,maxnp)
      common / dopt9 / dum1(maxnp),dum2(maxnp),dum3(maxnp),dum4(maxnp),
     *                 dum5(maxnp,maxnp),dum6(maxnp,maxnp),dum7(maxnp)
c
c     initialize size of blank common
c
      common calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns),
     *       cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns),
     *       omegax(maxnt,maxnt)
c
c
      write (*,1) maxalt,maxper,maxnt,maxnt1,
     *            maxnb,maxns,maxnp,maxnpp,
     *            maxbuf
1     format(
     *//' ============================================================'
     * /'       **** multinomial multiperiod probit model ****'
ccc     * /'    version 2.1  [sep 01, 1990]  (c) axel boersch-supan'
     * /'    version 2.1 - sep 01, 1990   (c) axel boersch-supan'
     * /' ============================================================'
     * /'       maxalt=',i3,' maxper=',i3,' maxnt=',i3,' maxnt1=',i3
     * /'       maxnb =',i3,' maxns =',i3,' maxnp=',i3,' maxnpp=',i3
     * /'                                         maxbuf=',i6
     * /' ============================================================'
     * /'      *** initialization, parameter and data input ***'
     * /' ============================================================')
      call inparm(parm,np,*998)
      call indata(*998)
c
c     ******************************
c     *   omega derivative check   *
c     ******************************
c
100   write (*,101)
101   format(
     *//' ==================================================='
     * /' multinomial multiperiod probit model:'
     * /'                              omega derivative check'
     * /' ===================================================')
      if (ipause.lt.0) then
         eps=0.000001
         if (mchk0.eq.0) goto 200
      else
         write (*,*) 'enter eps [1-9: 3=0.001 etc.; 0=no check] ?'
         read (*,'(i1)') mchk0
         if (mchk0.eq.0) goto 200
         if (mchk0.gt.0) eps=0.1**mchk0
      end if
c
      call dcomeg(parm,np,eps)
c
c     ****************************************************
c     *   initial probit/logit estimation (clark approx) *
c     ****************************************************
c
200   if (ipause.lt.0.and.method.eq.0) goto 300
201   write (*,202)
202   format(
     *//' ==============================================='
     * /' multinomial multiperiod probit model:'
     * /'                              initial estimation'
     * /' ===============================================')
      if (ipause.ge.0) then
         write (*,*) 
     *   ' enter initial estimation method [0=none, 1=logit, 2=clark] ?'
         read (*,'(i1)') method
      end if
      if (method.le.0.or.method.ge.4) goto 300
      if (method.eq.1) write (*,203)
      if (method.eq.2) write (*,204)
203   format(
     *//' ========================'
     * /' scaled multinomial logit'
     * /' ========================')
204   format(
     *//' =========================='
     * /' clark approximation probit'
     * /' ==========================')
      if (ispool.ne.0) write (kanal0,202)
      if (ispool.ne.0.and.method.eq.1) write (kanal0,203)
      if (ispool.ne.0.and.method.eq.2) write (kanal0,204)
c
c     **** derivative check
c
      ival=0
      idiff2=2
      if (ipause.ge.0) then
         write (*,*) ' derivative check of initial estimator',
     *                                    ' [0=none, 1=func, 2=check] ?'
         read (*,'(i1)') mchk1
      end if
      if (mchk1.ne.0) then
         if (ispool.ne.0) write (kanal0,'(//a)')
     *                 ' derivative check of initial estimator'
         call func(parm,np,fval,*998)
         write (*,'(a,f15.7)') ' fval=',fval
         if (ispool.ne.0) write (kanal0,'(a,f15.7)') ' fval=',fval
         if (mchk1.eq.2) then
            call deriv(parm,np,fval,grad,hess,idiff2,*998)
            write (*,'(a,30f10.4)') ' deriv=',(grad(i),i=1,np)
            if (ispool.ne.0) write (kanal0,'(a,30f10.4)')
     *                              ' deriv=',(grad(i),i=1,np)
            call fdiff(parm,np,fval,grad,*998)
            write (*,'(a,30f10.4)') ' fdiff=',(grad(i),i=1,np)
            if (ispool.ne.0) write (kanal0,'(a,30f10.4)')
     *                              ' fdiff=',(grad(i),i=1,np)
         end if
         call snooze(ipause)
      end if
c
c     **** maximize likelihood
c
      mx=1
      ialg=ialg1
      iter=iter1
      acc =accu1
      call opt(parm,grad,hess,hinv,np,
     *         xlf,fpnorm,grdir,ier,
     *         mx,iter,acc,ialg)
      call smlout(parm,hinv,np,
     *            xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul)
      call snooze(ipause)
      k0=0
      call promeg(parm,np,k0)
      if (ispool.ne.0) call promeg(parm,np,kanal0)
      call snooze(ipause)
c
c     *************************************************
c     *   method of simulated likelihood estimation   *
c     *************************************************
c
300   if (ipause.lt.0.and.nsimul.eq.0) goto 400
      write (*,301)
301   format(
     *//' ============================================================'
     * /' multinomial multiperiod probit model:'
     * /'                       smooth simulated likelihood estimation'
     * /' ============================================================')
      if (ipause.ge.0) then
         write (*,*) ' enter number of replications [0=quit] ?'
cccvaherror!!!         read (*,'(i1)') nsimul
         read (*,'(i4)') nsimul
         print *,nsimul
         call snooze(1)
      end if
      if (nsimul.eq.0) goto 400
      if (ispool.ne.0) write (kanal0,301)
c
      method=0
c
c     **** derivative check
c
      ival=0
      idiff2=2
      if (ipause.ge.0) then
          write (*,*) ' derivative check of sml estimator',
     *                                    ' [0=none, 1=func, 2=check] ?'
         read (*,'(i1)') mchk2
      end if
      if (mchk2.ne.0) then
         if (ispool.ne.0) write (kanal0,'(//a)')
     *                 ' derivative check of sml estimator'
         call func(parm,np,fval,*998)
         write (*,'(a,f15.7)') ' fval=',fval
         if (ispool.ne.0) write (kanal0,'(a,f15.7)') ' fval=',fval
         if (mchk2.eq.2) then
            call deriv(parm,np,fval,grad,hess,idiff2,*998)
            write (*,'(a,30f10.4)') ' deriv=',(grad(i),i=1,np)
            if (ispool.ne.0) write (kanal0,'(a,30f10.4)')
     *                              ' deriv=',(grad(i),i=1,np)
            call fdiff(parm,np,fval,grad,*998)
            write (*,'(a,30f10.4)') ' fdiff=',(grad(i),i=1,np)
            if (ispool.ne.0) write (kanal0,'(a,30f10.4)')
     *                              ' fdiff=',(grad(i),i=1,np)
         end if
         call snooze(ipause)
      end if
c
c     **** maximize likelihood
c
      mx=1
      ialg=ialg2
      iter=iter2
      acc =accu2
      call opt(parm,grad,hess,hinv,np,
     *         xlf,fpnorm,grdir,ier,
     *         mx,iter,acc,ialg)
      call smlout(parm,hinv,np,
     *            xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul)
      call snooze(ipause)
      k0=0
      call promeg(parm,np,k0)
      if (ispool.ne.0) call promeg(parm,np,kanal0)
c
c     *****************************
c     *   covariance correction   *
c     *****************************
c
400   if (ipause.lt.0.and.nwesml.eq.0) goto 900
      write (*,401)
401   format(
     *//' ============================================================'
     * /' multinomial multiperiod probit model:'
     * /'                       robust (wesml) covariance computation'
     * /' ============================================================')
      if (ipause.ge.0) then
         write (*,*) 'robust covariances ',
     *                              '[0=none, 1=old hess, 2=new hess] ?'
         read (*,'(i1)') nwesml
      end if
      if (nwesml.eq.0) goto 900
c
      if (ispool.ne.0) write (kanal0,401)
      call cwesml(parm,grad,hess,hinv,np,nwesml,*900)
      if (nwesml.eq.2) then
         ier1=100+ier
         call smlout(parm,hess,np,
     *               xlf,ier1,iter,ival,ialg,fpnorm,grdir,nsimul)
      end if
      ier2=200+ier
      call smlout(parm,hinv,np,
     *            xlf,ier2,iter,ival,ialg,fpnorm,grdir,nsimul)
      call snooze(ipause)
c
900   stop '*** sml finished ***'
998   stop '*** sml aborted ***'
      end
c
c
c
      subroutine func(parm,np,xlf,*)
c-----------------------------------------------------------------------
c     accumulates criterion function over individuals
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      dimension parm(np)
      integer*4 nn1,iseed,iseed1
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / dopt1 / parm2(maxnp),grad2(maxnp),hess2(maxnp,maxnp),
     *                 hinv2(maxnp,maxnp)
      common / ddiff / idiff2
      common / dival / ival
      common / doutp / kanal0,ispool,ipause
      common / dseed / iseed
      common / dsimu / nsimul
c
ccccccc    write (*,'(//a,i3)') '$$$ func: np,parm=',np
ccccccc    write (*,'(9f8.5)') (parm(i),i=1,np)
c
      if (method.eq.0.and.nsimul.eq.0) goto 98
c
      kobs=0
      nerr=0
      lderiv=.false.
      call omega1(parm,np,*97)
c     ---accumulate likelihood over nindiv
      nn1=1
      iseed1=iseed
      xlf=0.d0
      do 100 iobs=1,nindiv
         kobs=kobs+1
         call pull(nn1)
         call cont(xlf2,grad2,lderiv,iseed1,*95)
         xlf=xlf+xlf2
100      continue
c
c     error handling, exits
c
90    ival=ival+1
      return
95    write (*,*) ' func: failure in obs',kobs
      if (ispool.ne.0) write (kanal0,*) ' func: failure in obs',kobs
      return 1
97    write (*,*) ' func: failure in omega'
      if (ispool.ne.0) write (kanal0,*) ' func: failure in omega'
      return 1
98    write (*,*) ' func: misspecified nsimul=0'
      if (ispool.ne.0) write (kanal0,*) ' func: misspecified nsimul=0'
99    return 1
      end
c
      subroutine deriv(parm,np,xlf,grad,hess,idiff2,*)
c-----------------------------------------------------------------------
c
c     evaluates function, gradient (and hessian) of the
c     loglikelihood-function
c
c     idiff2=4:  bhhh approximation of hessian is cumulated
c                [hess = -grad*grad' = - inv(cov)]
c     otherwise: only function and gradient is computed
c                [hess will not be changed]
c
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      dimension parm(np),grad(np),grad2(maxnp),hess(np,np)
      integer*4 nn1,iseed,iseed1
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / dival / ival
      common / doutp / kanal0,ispool,ipause
      common / dseed / iseed
c
ccc      xlf=0.d0
      do 5 i=1,np
5        grad(i)=0.d0
      if (idiff2.eq.4) then
         do 6 i=1,np
         do 6 j=1,np
6           hess(i,j)=0.d0
      end if
      kobs=0
      lderiv=.true.
c     ---transfer parm on parmp and omega
      call omega2(parm,np,*95)
c     ---accumulate derivatives over nindiv
      nn1=1
      iseed1=iseed
      do 100 iobs=1,nindiv
         kobs=kobs+1
         call pull(nn1)
         call cont(xlf2,grad2,lderiv,iseed1,*95)
ccc         xlf=xlf+xlf2
         do 110 k=1,np
110         grad(k) = grad(k) + grad2(k)
         if (idiff2.eq.4) then
            do 150 j=1,np
            do 150 k=1,np
150            hess(j,k) = hess(j,k) - grad2(j)*grad2(k)
         end if
100   continue
c
c     error handling, exits
c
90    ival=ival+2
      return
95    write (*,*) ' failure in gradient evaluation, observation ',kobs
      if (ispool.ne.0) write (kanal0,*)
     *            ' failure in gradient evaluation, observation ',kobs
      return 1
      end
c
c
      subroutine cont(xlfi,grad,lderiv,iseed1,*)
c-----------------------------------------------------------------------
c
c     likelihood contribution xlfi for observation iobs
c
c     ** version for (nalt-1)*nper stochastic utility differences **
c     ** and first derivatives in grad
c
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      real*8 grad(maxnp)
      integer*4 iseed1
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / dxvar / x(maxnt,maxnpp),ic(maxper),wt
      common / dzbar / zquer(maxnt1),covz(maxnt1,maxnt1),
     *                 dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns)
c
c---------------------------------------------------
c    i  definition der zquer,
c       dzquer traegt ableitungen nach freien betas
c---------------------------------------------------
c
      it =0
      it1=0
      do 10 ip=1,nper
ccc                     write (*,*) '$$$ iper,nper=',ip,nper
         ialt=ic(ip)
         ita=it+ialt
         do 11 i=1,nalt
            it=it+1
            if (i.eq.ialt) goto 11
            it1=it1+1
ccc                           write (*,*) '$$$ i,ialt=',i,ialt
            sum=0
            do 12 k=1,nxy
ccc                           write (*,*) '$$$ it,it1,italt=',it,it1,ita
               xdiff=x(it,k)-x(ita,k)
               sum=sum+xdiff*parmp(k)
               i1=loc(k)
               if (i1.gt.0) dzquer(it1,i1)=xdiff
12          continue
            zquer(it1)=sum
11       continue
10    continue
c
ccc      write (*,'(//a)') '$$$ zquer='
ccc      write (*,'(9f8.5)') (zquer(i),i=1,nt1)
ccc      do 13 k=1,nb
ccc         write (*,'(//a)') '$$$ dzquer='
ccc13       write (*,'(9f8.5)') (dzquer(i,k),i=1,nt1)
c
c----------------------------------------------------
c     ii berechnung der covarianzmatrix von z,
c        dcovz traegt ableitungen nach freien sigmas
c----------------------------------------------------
c
      do 50 ip=1,nper
         ialt=ic(ip)
         ioff1=(ip-1)*nalt
         ioff2=(ip-1)*nalt1
         iii=ioff1+ialt
         do 50 jp=1,nper
            jalt=ic(jp)
            joff1=(jp-1)*nalt
            joff2=(jp-1)*nalt1
            jjj=joff1+jalt
c
            omegij=omega(iii,jjj)
c
            ki=1
            do 20 i=1,nalt
               if (i.ne.ialt) then
                  kj=1
                  do 25 j=1,nalt
c
                     if (j.ne.ialt) then
                        covz(ioff2+ki,joff2+kj)=
     *                            omega(ioff1+i,joff1+j)
     *                           -omega(ioff1+i,jjj    )
     *                           -omega(iii    ,joff1+j)
     *                           +omegij
c
                        if (lderiv) then
                           do 26 is=1,ns
                              dcovz(ioff2+ki,joff2+kj,is)=
     *                            domega(ioff1+i,joff1+j,is)
     *                           -domega(ioff1+i,jjj    ,is)
     *                           -domega(iii    ,joff1+j,is)
     *                           +domega(iii    ,jjj    ,is)
26                         continue
                        end if
c
                        kj=kj+1
                     end if
c
25                continue
                  ki=ki+1
               end if
20          continue
c
50    continue
c
ccc      write (*,'(//a)') '$$$ omega='
ccc      do 58 i=1,nt
ccc58       write (*,'(9f8.5)') (omega(i,j),j=1,nt)
ccc      write (*,'(//a)') '$$$ covz='
ccc      do 59 i=1,nt1
ccc59       write (*,'(9f8.5)') (covz(i,j),j=1,nt1)
c
c------------------------------------------------------
c     iii  probability evaluation
c-------------------------------------------------------
c
900   if (method.eq.1) then
       call logit2(zquer,dzquer,prob,grad,lderiv,*991)
      else if (method.eq.2) then
       call clark2(zquer,dzquer,covz,dcovz,prob,grad,lderiv,*991)
      else
       call gewek2(zquer,dzquer,covz,dcovz,prob,grad,lderiv,iseed1,*991)
      end if
c
      if (prob.lt.1.d-300) prob=1.d-300
cccc      write (*,*) '$$$ average prob=',prob
c
      xlfi=wt*dlog(prob)
      if (lderiv) then
      do 910 k=1,nb+ns
910      grad(k)=wt*grad(k)/prob
      end if
c
      return
991   return 1
      end
c
c
c
c======================================================================
c     memory resident utilities
c======================================================================
c
c
c
      subroutine snooze(ipause)
c----------------------------------------------------------
c     gives you a little break to view output
c
c     ipause<0 skips
c     ipause in seconds<60. requires system clock call.
c     ipause=60 pauses
c----------------------------------------------------------
      implicit integer (i-n)
      integer*4 isec0,isec1,isec2
      integer*2 i0,i1,i2,i3,i4
      if (ipause.le.0) goto 90
      if (ipause.ge.60) goto 10
ccc      write (*,'(a,i2,a\)') ' ...pausing for ',ipause,' seconds'
      write (*,'(a,i2,a)') ' ...pausing for ',ipause,' seconds'
      i0=0
      call gettim(i1,i2,i0,i4)
      isec0=i1*3600+i2*60+i0
      isec2=isec0
      isec=0
      i3=0
      do 1000 i=1,32000
         call gettim(i1,i2,i3,i4)
         isec1=i1*3600+i2*60+i3
         isec=isec1-isec0
         if (isec1.ne.isec2) then
ccc            write (*,'(a1\)') '.'
            write (*,'(a1)') '.'
            isec2=isec1
         end if
         if (isec.gt.ipause) goto 90
1000  continue
      goto 90
10    write (*,'(/'' hit enter to continue.....''/)')
      read (*,'(a1)')
90    write (*,'(a1)') ' '
      return
      end
c
      subroutine tstamp(tdate)
c----------------------------------------------------------
c     generates string*21 tdate with time and date
c     e.g.: 31 dec 1989, 12:00:00
c----------------------------------------------------------
      implicit integer (i-n)
      integer*2 iyear,imon,iday,ihour,imin,isec,ihsec
      character*21 tdate
      character*3 month(12)
      character*2 s1,s2,s3,s4,s5
      data month/'jan','feb','mar','apr','may','jun',
     *           'jul','aug','sep','oct','nov','dec'/
      call getdat(iyear,imon,iday)
      isec=0
      call gettim(ihour,imin,isec,ihsec)
      call int2ch(iday,s1)
      iyear=iyear-1900
      call int2ch(iyear,s2)
      call int2ch(ihour,s3)
      call int2ch(imin,s4)
      call int2ch(isec,s5)
      tdate=s1//' '//month(imon)//' 19'//s2//', '//s3//':'//s4//':'//s5
      return
      end
c
      subroutine int2ch(int,ch2)
c----------------------------------------------------------
c     converts 2-digit integer to 2-byte character
c----------------------------------------------------------
      integer*2 int,i1,i2
      character*1 c1,c2
      character*2 ch2
      i1=int/10
      i2=int-i1*10
      if (i1.le.9.and.i2.le.9.and.i1.ge.0.and.i2.ge.0) then
         c1=char(i1+48)
         c2=char(i2+48)
         ch2=c1//c2
      else
         ch2='??'
      end if
      return
      end
c
      subroutine openf(kanal,text,defaul,mode)
c---------------------------------------------------------------------
c     conditional and interactive opening of files
c     mode:  1 = old, formatted     2 = old, unformatted
c            3 = new, formatted     4 = new, unformatted
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      character prompt*1,defaul*40,text*40,infile*40
      logical exist,fmt
      common / doutp / kanal0,ispool,ipause
c
10    write (*,11) text,defaul
11    format(//' enter filename for ',a40
     *        /' [default name is ',a40,']'
ccc     *        /' >'\)
     *        /' >')
      if (ipause.ge.0) then
         read (*,'(a40)') infile
         if (infile.eq.' ') infile=defaul
         if (infile.eq.' ') goto 10
      else
         infile=defaul
      end if
      write (*,12) infile,text
      if (ispool.ne.0.and.kanal.ne.kanal0) write (kanal0,12) infile,text
12    format(/' opening ',a40
     *       /'     for ',a40)
      inquire (file=infile,exist=exist)
c
      if (mode.eq.1.or.mode.eq.3) then
         fmt=.true.
      else
         fmt=.false.
      end if
c
      if (mode.le.2) then
         if (exist) then
            if (fmt) then
               open (kanal,file=infile,status='old',form='formatted')
            else
               open (kanal,file=infile,status='old',form='unformatted')
            end if
         else
            write (*,*) '*** error, does not exist'
            if (ipause.lt.0) stop '*** hlogit aborted ***'
            goto 10
         end if
      else
         if (exist) then
            write (*,'('' *** warning: overwrite ? (y/n) >'')')
            if (ipause.ge.0) then
               read  (*,'(a1)') prompt
               if (prompt.eq.'n'.or.prompt.eq.'n') goto 10
            end if
            if (fmt) then
               open (kanal,file=infile,status='old',form='formatted')
            else
               open (kanal,file=infile,status='old',form='unformatted')
            end if
         else
            if (fmt) then
               open (kanal,file=infile,status='new',form='formatted')
            else
               open (kanal,file=infile,status='new',form='unformatted')
            end if
         end if
      end if
      return
      end
c
      subroutine matp(a,n,nadim,k,b,nbdim,l,c,ncdim)
c----------------------------------------------------------------------
c     c(n,l) = a(n,k) * b(k,l)
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension a(nadim,k),b(nbdim,l),c(ncdim,l)
c
      do 50 i=1,n
      do 50 j=1,l
         sum=0.d0
         do 51 m=1,k
51          sum=sum+a(i,m)*b(m,j)
50       c(i,j)=sum
      return
      end
c
      real*8 function dotp(x,np,y)
c----------------------------------------------------------------------
c     dotp=x'y
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension x(np),y(np)
      dotp=0.
      do 10 i=1,np
 10      dotp=dotp+x(i)*y(i)
      return
      end
c
      subroutine move(x1,np,x2)
c----------------------------------------------------------------------
c     vector x2=x1
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension x1(np),x2(np)
      do 10 i=1,np
10       x2(i)=x1(i)
      return
      end
c
      subroutine movem(a1,a2,np,maxnp)
c----------------------------------------------------------------------
c     matrix a2=a1
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension a1(maxnp,np),a2(maxnp,np)
      do 10 i=1,np
      do 10 j=1,np
10       a2(i,j)=a1(i,j)
      return
      end
c
      subroutine step(parm,np,size,direc,pnew)
c----------------------------------------------------------------------
c     pnew=parm+size*direc
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension parm(np),direc(np),pnew(np)
      do 10 i=1,np
10       pnew(i)=parm(i)+size*direc(i)
      return
      end
c
      subroutine gaussj(a,b,n,npmax,*)
c---------------------------------------------------------------------
c     solves ax=b by gauss-jordan elimination.  full pivoting.
c     on return: a carries inv(a), b carries x
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      parameter (nmax=50)
      dimension a(npmax,n),b(n),ipiv(nmax),indxr(nmax),indxc(nmax)
      do 11 j=1,n
         ipiv(j)=0
11    continue
      do 22 i=1,n
         big=0.0
         do 13 j=1,n
            if (ipiv(j).ne.1) then
               do 12 k=1,n
                  if (ipiv(k).eq.0) then
                     if (abs(a(j,k)).ge.big) then
                        big=abs(a(j,k))
                        irow=j
                        icol=k
                     endif
                  else if (ipiv(k).gt.1) then
                     return 1
                  endif
12             continue
            endif
13       continue
         ipiv(icol)=ipiv(icol)+1
         if (irow.ne.icol) then
            do 14 l=1,n
               dum=a(irow,l)
               a(irow,l)=a(icol,l)
               a(icol,l)=dum
14          continue
            dum=b(irow)
            b(irow)=b(icol)
            b(icol)=dum
         endif
         indxr(i)=irow
         indxc(i)=icol
         if (a(icol,icol).eq.0.0) return 1
         pivinv=1./a(icol,icol)
         a(icol,icol)=1.0
         do 16 l=1,n
            a(icol,l)=a(icol,l)*pivinv
16       continue
         b(icol)=b(icol)*pivinv
         do 21 ll=1,n
            if (ll.ne.icol) then
               dum=a(ll,icol)
               a(ll,icol)=0.0
               do 18 l=1,n
                  a(ll,l)=a(ll,l)-a(icol,l)*dum
18             continue
               b(ll)=b(ll)-b(icol)*dum
            endif
21       continue
22    continue
      do 24 l=n,1,-1
         if (indxr(l).ne.indxc(l)) then
            do 23 k=1,n
               dum=a(k,indxr(l))
               a(k,indxr(l))=a(k,indxc(l))
               a(k,indxc(l))=dum
23          continue
         endif
24    continue
      return
      end
cc
      subroutine trunc2(urand,rr,eps,prob,deps,dprob)
c-----------------------------------------------------------------------
c     eps: right truncated normal at rr for uniform random variate urand
c     prob: prob(-inf < standard normal variate < rr)
c
c     deps,dprob: first derivatives of eps and prob with respect to rr
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      prob=cumphi(rr)
      dprob=denphi(rr)
      p=prob*urand
      dp=dprob*urand
      call phinv2(p,eps,deps)
      deps=deps*dp
      if (eps.le.-8.25) then
cover    write (*,'(a,f12.5,d25.17)')
cover*         ' overflow in cuminv: rr,h(rr)=',rr,prob
         if (eps.gt.rr) eps=rr
      end if
      return
      end
c
      subroutine trunc(urand,rr,eps,prob)
c-----------------------------------------------------------------------
c     eps: right truncated normal at rr for uniform random variate urand
c     prob: prob(-inf < standard normal variate < rr)
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      prob=cumphi(rr)
      p=prob*urand
      call phinv(p,eps)
      if (eps.le.-8.25) then
cover    write (*,'(a,f12.5,d25.17)')
cover*         ' overflow in cuminv: rr,h(rr)=',rr,prob
         if (eps.gt.rr) eps=rr
      end if
      return
      end
c
      subroutine phinv2(prob,arg,darg)
c---------------------------------------------------------------------
c
c     inverse normal c.d.f.: prob=cumphi(cuminv(prob))
c                            arg =cuminv(cumphi(arg ))
c
c     (1) formula 26.2.23 of abramovitz and stegun p.933 10th edition.
c         (accuracy<4.5e-4)
c     (2) inverse interpolation to 3rd order (a&s example 5 p.954)
c
c     note: "darg" carries first derivative of "arg" w.r.t. "prob".
c
c     modified from vah 24sep84/14apr85.
c     axel boersch-supan (c) version 25. august 1990
c
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      data c0,c1,c2/2.515517d0,.802853d0,.010328d0/,
     *     d1,d2,d3/1.432788d0,.189269d0,.001308d0/
c
c     crude ways of dealing with prob=0 and prob=1
c
      if (prob.lt.1.d-37) then
         arg=-8.25d0
         darg=2.d15
         return
      end if
      p1=1.d0-prob
      if (p1.lt.1.d-37) then
         arg=8.25d0
         darg=2.d15
      return
      end if
c
c     flip if prob>0.5
c
      p=prob
      iflip=0
      if (prob.gt.0.5d0) then
         p=p1
         iflip=1
      end if
c
c     p<0.5 as required by inversion:
c
      t=dsqrt(dlog(1.d0/(p**2)))
      dt=-1.d0/(t*p)
      top =c0+c1*t+c2*t**2
      dtop=   c1  +2*c2*t
      bot =1+d1*t+d2*t**2+d3*t**3
      dbot=  d1  +2*d2*t +3*d3*t**2
      x0=t-top/bot
      dx0=(1-dtop/bot+(dbot/bot)*(top/bot))*dt
c
c     up to now x0<0:
c
      x0=-x0
      dx0=-dx0
c
c     now x0>0 as required by inverse interpolation:
c
      denx=denphi(x0)
      cumx=cumphi(x0)
      ddenx=-denx*x0*dx0
      dcumx= denx*dx0
      t=(p-cumx)/denx
      dt=(1-dcumx)/denx-t*ddenx/denx
c
      ph1  = x0 + t + 0.5* x0*t**2
      dph1 =dx0 +dt + 0.5*dx0*t**2 + x0*t*dt
      ph2  = ((1+2*x0**2)*t**3)/6.d0
      dph2 = (   4*x0*dx0*t**3)/6.d0 +((1+2*x0**2)*t**2)*dt/2.d0
      arg=ph1+ph2
      darg=dph1+dph2
c
c     flip back if prob>0.5
c
      if (iflip.eq.1) arg=-arg
      return
      end
c
      subroutine phinv(prob,arg)
c---------------------------------------------------------------------
c
c     inverse normal c.d.f.: prob=cumphi(cuminv(prob))
c                            arg =cuminv(cumphi(arg ))
c
c     (1) formula 26.2.23 of abramovitz and stegun p.933 10th edition.
c         (accuracy<4.5e-4)
c     (2) inverse interpolation to 3rd order (a&s example 5 p.954)
c
c     modified from vah 24sep84/14apr85.
c     axel boersch-supan (c) version 25. august 1990
c
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      data c0,c1,c2/2.515517d0,.802853d0,.010328d0/,
     *     d1,d2,d3/1.432788d0,.189269d0,.001308d0/
c
c     crude ways of dealing with prob=0 and prob=1
c
      if (prob.lt.1.d-37) then
         arg=-8.25d0
         darg=2.d15
         return
      end if
      p1=1.d0-prob
      if (p1.lt.1.d-37) then
         arg=8.25d0
         darg=2.d15
      return
      end if
c
c     flip if prob>0.5
c
      p=prob
      iflip=0
      if (prob.gt.0.5d0) then
         p=p1
         iflip=1
      end if
c
c     p<0.5 as required by inversion:
c
      t=dsqrt(dlog(1.d0/(p**2)))
      top =c0+c1*t+c2*t**2
      bot =1+d1*t+d2*t**2+d3*t**3
      x0=t-top/bot
c
c     up to now x0<0:
c
      x0=-x0
c
c     now x0>0 as required by inverse interpolation:
c
      denx=denphi(x0)
      cumx=cumphi(x0)
      t=(p-cumx)/denx
c
      ph1  = x0 + t + 0.5* x0*t**2
      ph2  = ((1+2*x0**2)*t**3)/6.d0
      arg=ph1+ph2
c
c     flip back if prob>0.5
c
      if (iflip.eq.1) arg=-arg
      return
      end
c
c
c
c
      real*8 function denphi(x)
      implicit real*8 (a-h,o-z)
c     1/sqrt(2*pi)
      data constp/0.39894228040143267793d00/
      t=-0.5*x*x
      if (t.lt.-707.0) then
cover    write (*,*) 'denphi(x) underflow. x=',x
         denphi=0
      else
         denphi=dexp(t)*constp
      end if
      return
      end
c
      real*8 function cumphi(x)
      implicit real*8 (a-h,o-z)
      data sqrt2/1.414213562d0/
      if (x.ge.0) then
         z=x/sqrt2
         cumphi=1.0-0.5*erfc(z)
      else
         z=-x/sqrt2
         cumphi=0.5*erfc(z)
      end if
      return
      end
c
      real*8 function erfc(x)
      implicit real*8 (a-h,o-z)
      data a0/-1.26551223d0/,
     *     a1/ 1.00002368d0/,
     *     a2/ 0.37409196d0/,
     *     a3/ 0.09678418d0/,
     *     a4/-0.18628806d0/,
     *     a5/ 0.27886807d0/,
     *     a6/-1.13520398d0/,
     *     a7/ 1.48851587d0/,
     *     a8/-0.82215223d0/,
     *     a9/ 0.17087277d0/
      z=dabs(x)
      t=1.0/(1.0+0.5*z)
      p9=a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9))))))))
      zzz=-z*z+p9
      if (zzz.lt.-707.0) then
cover    write (*,*) 'erfc(x) underflow. x=',x
         zzz=0
      else
         zzz=t*exp(zzz)
      end if
      if (x.lt.0) zzz=2.0-zzz
      erfc=zzz
      return
      end
c
      real*8 function randu(iseed)
c-----------------------------------------------------------------------
c     generates randu pseudo uniform [0,1] random number generator.
c     requires a seed (iseed) which is updated and returned.
c     [lewis,goodman & miller,ibm systems journal,v.8#2 1969,pp.136-145]
c-----------------------------------------------------------------------
      implicit real*8 (a-z)
      integer*4 iseed
      data xb,xm/16807.d0,2147483647.d0/
      seed=iseed
      xseed=dmod(xb*seed,xm)
      randu=xseed/xm
      iseed=xseed
      return
      end
c
c
      subroutine chol2(a,da,c,dc,n,maxn,ns,maxns,*)
c----------------------------------------------------------------------
c     cholesky decomposition of symmetric nxn matrix a (dimension maxn)
c     results in lower triangular matrix c with a=c*c'
c
c     *** version with first partials: **
c     dc carries first partials of c  w.r.t.  ns deep parameters in a,
c     if da carries first partials of a  w.r.t.  deep parameters.
c
c     modified from stoer/bulirsch
c     (c) axel boersch-supan     version: 25. august 1990
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension a(maxn,maxn),c(maxn,maxn),
     *          da(maxn,maxn,maxns),dc(maxn,maxn,maxns),
     *          dx(10)
      notpos=0
      do 10 i=1,n
      do 10 j=i,n
         x=a(i,j)
         do 101 is=1,ns
101         dx(is)=da(i,j,is)
         do 20 kk=1,i-1
            k=i-kk
            x=x-a(j,k)*a(i,k)
            do 201 is=1,ns
201            dx(is)=dx(is)-da(j,k,is)*a(i,k)-a(j,k)*da(i,k,is)
20       continue
         if (i.eq.j) then
            if (x.le.0) goto 99
            sqrtx=dsqrt(x)
            c(i,i)=1.0/sqrtx
            do 211 is=1,ns
211            dc(i,i,is)=-0.5*dx(is)/(x*sqrtx)
         else
            a(j,i)=x*c(i,i)
            do 221 is=1,ns
221            da(j,i,is)=dx(is)*c(i,i)+x*dc(i,i,is)
         end if
10    continue
30    do 31 i=1,n
         do 32 j=i,n
            if (i.eq.j) goto 32
               c(j,i)=a(j,i)
               c(i,j)=0.0
               a(j,i)=a(i,j)
               do 301 is=1,ns
                  dc(j,i,is)=da(j,i,is)
                  dc(i,j,is)=0.0
301               da(j,i,is)=da(i,j,is)
32       continue
         dd=c(i,i)
         c(i,i)=1.0/dd
         dd=dd**2
         do 321 is=1,ns
321         dc(i,i,is)=-dc(i,i,is)/dd
31    continue
      if (notpos.eq.1) return 1
      goto 90
99    write (*,*) 'cholesky: matrix not positive definite'
      notpos=1
      goto 30
90    return
      end
c
      subroutine chol(a,c,n,maxn,*)
c----------------------------------------------------------------------
c     cholesky decomposition of symmetric nxn matrix a (dimension maxn)
c     results in lower triangular matrix c with a=c*c'
c
c     modified from stoer/bulirsch
c     (c) axel boersch-supan     version: 25. august 1990
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension a(maxn,maxn),c(maxn,maxn)
      notpos=0
      do 10 i=1,n
      do 10 j=i,n
         x=a(i,j)
         do 20 kk=1,i-1
            k=i-kk
            x=x-a(j,k)*a(i,k)
20       continue
         if (i.eq.j) then
            if (x.le.0) goto 99
            sqrtx=dsqrt(x)
            c(i,i)=1.0/sqrtx
         else
            a(j,i)=x*c(i,i)
         end if
10    continue
30    do 31 i=1,n
         do 32 j=i,n
            if (i.eq.j) goto 32
               c(j,i)=a(j,i)
               c(i,j)=0.0
               a(j,i)=a(i,j)
32       continue
         dd=c(i,i)
         c(i,i)=1.0/dd
31    continue
      if (notpos.eq.1) return 1
      goto 90
99    write (*,*) 'cholesky: matrix not positive definite'
      notpos=1
      goto 30
90    return
      end
c
      subroutine pull(nn1)
c-----------------------------------------------------------------------
c     pulls data for one observation from internal worksheet.
c     (nda items beginning with offset nn1)
c
c     this routine in combination with sr trans is designed to be the
c     only access to the internal worksheet.
c     the nesting (pull calls trans) is necessary for the data(*) call
c     in order to isolate routines that must be compiled with the $large
c     metacommand (under ms-fortran).
c
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      real*4  data(maxbuf)
      integer*4 nn1
      common / ddata / data
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      call trans(data(nn1))
      nn1=nn1+ndp
      return
      end
c
ccc$large
      subroutine trans(data)
c-----------------------------------------------------------------------
c     transfers relevant piece of data(*) on matrix x in common dtrans.
c     if maxbuf is larger than one memory segment, this (and only this)
c     routine must be compiled with the $large metacommand (ms-fortran)
c-----------------------------------------------------------------------
c     data set-up:  x-vars (nx), choice (1), weight(nwt), y-vars (ny)
c     ndd=nx*nalt
c     ndw=nx*nalt+1+nwt
c     nda=nx*nalt+1+nwt+ny  = number of data items per observation
c     ndp=nda*nper          = number of data items per individual
c     nb=nx+nalt1*(ny+nc)   = number of parameters in beta/x
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c-----$large data----------------------include, if maxbuf>16000---------
      real*4  data(*)
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dxvar / x(maxnt,maxnpp),ic(maxper),wt
c
ccccc write (*,'(a,7f9.4,5(/8f9.4))') '$$$ data=',(data(i),i=1,ndp)
c
      mper1=0
      mper2=0
      do 1000 iper=1,nper
c
c        chosen alternative and weight
c
         ic(iper)=data(mper1+ndd+1)+.5
                ialt=ic(iper)
                if (ialt.lt.1.or.ialt.gt.nalt) then
                  write (*,*) 'error in cont: ialt=',data(mper1+ndd+1)
                  write (*,*) '  iper,ioff,ndd=',iper,mper1+ndd+1,ndd
                end if
         wt=1.0
         if (nwt.eq.1) wt=data(mper1+ndd+2)
c
c        transfer observation-relevant piece of data(*)
c
c        ---alternative-specific variables
         k1=0
         do 10 k=1,nx
         do 10 i=1,nalt
            k1=k1+1
10          x(mper2+i,k)=data(mper1+k1)
c        ---agent-specific variables
         k1=nx
         do 20 k=1,ny
            da=data(mper1+ndw+k)
            do 20 k2=1,nalt1
               k1=k1+1
               do 22 i=1,nalt
22                x(mper2+i,k1)=0.0
               x(mper2+k2,k1)=da
20       continue
c        ---alternative-specific constants
         do 30 k2=1,nalt1
            k1=k1+1
            do 32 i=1,nalt
32             x(mper2+i,k1)=0.0
            x(mper2+k2,k1)=1.0
30       continue
c
         mper1=mper1+nda
         mper2=mper2+nalt
1000     continue
c
cc      i=0
cc      do 50 i1=1,nper
cc         write (*,*) 'nper,iper,ic(iper) / x=',nper,i1,ic(i1)
cc         do 50 i2=1,nalt
cc            i=i+1
cc50          write (*,'(a,7f9.4/8f9.4)') '$$$ x=',(x(i,k),k=1,nxy)
c
      return
      end
c
c$include smlfunc.for
c$include smlinp.for
c$include smlutil.for
c$include smlopt.for
c$include wf.add
c---------------------------------------------------------------------
c.pl100/.rm72/
c     sml-func: likelihood function
c     =============================
c
c     copyright (c) axel boersch-supan,  version 10. march 1990
c
c     *** updated from mppmss.for:       1. july   1990
c     *** updated from mslsiml.for:      27. august 1990
c     *** merged from sml-siml/init.for: 1. september 1990
c
c     subroutines:  gewek2
c                   clark2
c                   logit2
c                   omega1
c                   omega2
c
c---------------------------------------------------------------------
c
      subroutine gewek2(zquer,dzquer,covz,dcovz,
     *                                      prob,dprob,lderiv,iseed1,*)
c-----------------------------------------------------------------------
c
c     geweke algorithm to generate unbiased nt-dim choice probability
c
c     u = n(0,covz)   s.t.  .mdnm/-l s u s zquer
c
c
c     axel boersch-supan  (c)  25. august 1990
c
c     ***  version with analytic first derivatives dprob 
c                       w.r.t. deep parms in covz
c
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      real*8 dprob(maxnp),
     *        zquer(maxnt1),covz(maxnt1,maxnt1),
     *        dzquer(maxnt1,maxnb),dcovz(maxnt1,maxnt1,maxns)

      integer*4 iseed1
      logical   lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dsimu / nsimul
      common / dfunc / cmat(maxnt1,maxnt1),e(maxnt1),
     *                 dcmat(maxnt1,maxnt1,maxns),de(maxnt1,maxnp),
     *                 dgprob(maxnp),drr(maxnp)
      common / doutp / kanal0,ispool,ipause
c
c     note: w-dimension = nt1 = (nalt-1)*nper
c
c     initialize
c
      prob=0
      np=nb+ns
      mxt1=maxnt1
c
c     *** get cholesky factor of covz (on lower triangle)
c
      if (lderiv) then
         do 10 ip=1,np
10          dprob(ip)=0
         mxs=maxns
         call chol2(covz,dcovz,cmat,dcmat,nt1,mxt1,ns,mxs,*990)
      else
         call chol(covz,cmat,nt1,mxt1,*990)
      end if
c
c     loop on nsimul
c
c     *** draw truncated normal eps and associated probability for each
c         alternative and period, store on e and probs
c
      do 100 isimul=1,nsimul
c
         j=1
         dd=cmat(1,1)
         rr=-zquer(1)/dd
         urand=randu(iseed1)
cccc         write (*,*) 'isimul,urand=',isimul,urand
c
         if (lderiv) then
            call trunc2(urand,rr,eps,zprob,deps,dzprob)
            e(1)=eps
            gprob=zprob
            do 110 ib=1,nb
110            drr(ib)=-dzquer(1,ib)/dd
            do 111 is=1,ns
               ip=nb+is
111            drr(ip)=-rr*dcmat(1,1,is)/dd
            do 115 ip=1,np
               de(1,ip)=deps*drr(ip)
115            dgprob(ip)=dzprob*drr(ip)
         else
            call trunc(urand,rr,eps,zprob)
            e(1)=eps
            gprob=zprob
         end if
c
         do 200 j=2,nt1
            dd=cmat(j,j)
            rr=0
            do 210 i=1,j-1
210            rr=rr+cmat(j,i)*e(i)
            rr=(-zquer(j)-rr)/dd
            urand=randu(iseed1)
c
            if (lderiv) then
               call trunc2(urand,rr,eps,zprob,deps,dzprob)
               e(j)=eps
               oldgpr=gprob
               if (zprob.lt.1.d-137) zprob=0.d0
               gprob=gprob*zprob
               do 220 ib=1,nb
                  drr(ib)=0
                  do 221 i=1,j-1
221                  drr(ib)=drr(ib)+cmat(j,i)*de(i,ib)
                  drr(ib)=(-dzquer(j,ib)-drr(ib))/dd
220            continue
               do 222 is=1,ns
                  ip=nb+is
                  drr(ip)=0
                  do 223 i=1,j-1
223                   drr(ip)=drr(ip)+dcmat(j,i,is)*e(i)
     *                               + cmat(j,i)   *de(i,ip)
                  drr(ip)=-drr(ip)/dd-rr*dcmat(j,j,is)/dd
222            continue
               do 225 ip=1,np
                  de(j,ip)=deps*drr(ip)
                  dgprob(ip)=dgprob(ip)*zprob+oldgpr*dzprob*drr(ip)
225            continue
            else
               call trunc(urand,rr,eps,zprob)
               e(j)=eps
               if (zprob.lt.1.d-137) zprob=0.d0
               gprob=gprob*zprob
            end if
c
200      continue
ccccccc              write (*,*) 'isimul,prob=',isimul,gprob
c
c        accumulate probability
c
         prob=prob+gprob
         if (lderiv) then
            do 120 ip=1,np
120            dprob(ip)=dprob(ip)+dgprob(ip)
         end if
c
c     end loop on nsimul
c
100   continue
c
c     average across nsimul
c
900   prob=prob/float(nsimul)
      if (lderiv) then
         do 901 ip=1,np
901         dprob(ip)=dprob(ip)/float(nsimul)
      end if
c
      return
990   write (*,*) '*** error in geweke: covz not pos.definite ***'
      if (ispool.ne.0) write (kanal0,*)
     *            '*** error in geweke: covz not pos.definite ***'
      return 1
      end
c
c
      subroutine omega1(parm,np,*)
c=======================================================================
c
c     creates covariance matrix omega (abbreviated from omega2)
c
c     (c) axel boersch-supan,  27. august 1990
c
c======================================================================
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c
      dimension parm(np)
c
      common calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns),
     *       cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns),
     *       omegax(maxnt,maxnt)
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / doutp / kanal0,ispool,ipause
c
c     (a) transfer updated parm vector on full parmp vector
c         and initialize matrices with identity/zero matrix
c
      do 100 i=1,np
100      parmp(loc1(i))=parm(i)
c
      do 110 i=1,nt
         do 111 j=1,nt
            omega(i,j)=0.0
111      continue
         omega(i,i)=1.0
110   continue
c
c     (b) correlation matrix of unobserved utility components
c         (cases 2,4)
c
2000  if (model.ne.2.and.model.ne.4) goto 3000
c
      do 210 ia=1,nalt
      do 210 ja=1,nalt
c
c        fetch standard deviation and correlation of unobserved
c        utility components and undo transformation
c
         ns1i=ns1+ia
         ns1j=ns1+ja
         ns2ij=ns2+maps(ia,ja)
         stni=parmp(ns1i)
         if (dabs(stni).gt.707.0) goto 9000
         stni=dexp(stni)
         stnj=parmp(ns1j)
         if (dabs(stnj).gt.707.0) goto 9000
         stnj=dexp(stnj)
         if (ia.eq.ja) then
            cnij=1.0
         else
            temp=parmp(ns2ij)
            if (dabs(temp).gt.707.0) goto 9000
            temp=dexp(temp)
            cnij=(temp-1)/(temp+1)
         end if
c
c        compute covariance of unobserved utilities
c
         cnij=stni*stnj*cnij
         cnunu(ia,ja)=cnij
c
210   continue
c
c     (c) correlation matrix of random effects
c         (cases 3,4,7,8)
c
3000  if (model.le.2.or.model.eq.5.or.model.eq.6) goto 4000
c
      do 310 ia=1,nalt
      do 310 ja=1,nalt
c
c        fetch standard deviation and correlation of random effects
c        and undo transformation
c
         ns3i=ns3+ia
         ns3j=ns3+ja
         ns4ij=ns4+maps(ia,ja)
         stai=dexp(parmp(ns3i))
         staj=dexp(parmp(ns3j))
         if (ia.eq.ja) then
            caij=1.0
         else
            temp=dexp(parmp(ns4ij))
            caij=(temp-1)/(temp+1)
         end if
c
c        compute covariance of unobserved utilities
c
         caij=stai*staj*caij
         calfa(ia,ja)=caij
c
310   continue
c
c     (d) autoregressive process with iid nu's
c         (cases 5,7)
c
4000  if (model.ne.5.and.model.ne.7) goto 5000
c
      do 400 ia=1,nalt
      do 400 ja=1,nalt
c
c        fetch rho and undo transformation
c
         ns5i=ns5+ia
         ns5j=ns5+ja
         temp=dexp(parmp(ns5i))
         rhoi=(temp-1)/(temp+1)
         temp=dexp(parmp(ns5j))
         rhoj=(temp-1)/(temp+1)
c
c        identity for nu's
c
         cnij=0.0
         if (ia.eq.ja) cnij=1.0
c
c        time dependent part
c
         do 420 it=1,nper
         do 420 jt=1,nper
            is=it-jt
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
            if (jt.lt.it.and.rhoi.ne.0) then
               rhoit=rhoi**is
               omega(ita,jta)=rhoit*cnij
            else if (it.lt.jt.and.rhoj.ne.0) then
               rhojt=rhoj**(-is)
               omega(ita,jta)=rhojt*cnij
            else if (it.eq.jt) then
               omega(ita,jta)=cnij
            end if
420      continue
400   continue
c
c     (e) autoregressive process with correlated nu's
c         (cases 6,8)
c
5000  if (model.ne.6.and.model.ne.8) goto 8000
c
      do 500 ia=1,nalt
      do 500 ja=1,nalt
c
c        fetch std.dev's and correlation of unobserved util's
c        and undo transformation
c
         ns1i=ns1+ia
         ns1j=ns1+ja
         ns2ij=ns2+maps(ia,ja)
         stni=dexp(parmp(ns1i))
         stnj=dexp(parmp(ns1j))
         if (ia.eq.ja) then
            cnij=1.0
         else
            temp=dexp(parmp(ns2ij))
            cnij=(temp-1)/(temp+1)
         end if
c
c        fetch rho and and undo transformations
c
         ns5i=ns5+ia
         ns5j=ns5+ja
         temp=dexp(parmp(ns5i))
         rhoi=(temp-1)/(temp+1)
         temp=dexp(parmp(ns5j))
         rhoj=(temp-1)/(temp+1)
         rho1ij=1.0-rhoi*rhoj
c
c        stationarity assumption
c
         snij=stni*stnj/rho1ij
         cnij=snij*cnij
c
c        time dependent part
c
         do 520 it=1,nper
         do 520 jt=1,nper
            is=it-jt
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
c
c           3 cases: it>jt (upper right blocks)
c                    it=jt (diagonal blocks)
c                    it<jt (lower left blocks)
c
            if (jt.lt.it.and.rhoi.ne.0) then
               rhos=rhoi**is
               omega(ita,jta)=rhos*cnij
            else if (it.lt.jt.and.rhoj.ne.0) then
               rhos=rhoj**(-is)
               omega(ita,jta)=rhos*cnij
            else if (it.eq.jt) then
               omega(ita,jta)=cnij
            end if
520      continue
500   continue
c
c     (g) add pieces together, according to case
c
8000  goto (8100,8200,8300,8400,8100,8100,8500,8500),model
c            1    2    3    4    5    6    7    8
c
c     model=1,5,6: done
c
8100  return
c
c     model=2: transfer cnunu block to all diagonal blocks
c
8200  do 8250 ip=1,nper
         ipt=(ip-1)*nalt
         do 8260 ia=1,nalt
         do 8260 ja=1,nalt
            id=ipt+ia
            jd=ipt+ja
            omega(id,jd)=cnunu(ia,ja)
8260     continue
8250  continue
      return
c
c     model=3: calfa on all blocks plus identity matrix on diag. blocks
c
c     transfer calfa block to all blocks
c
8300  do 8350 ip=1,nper
      do 8350 jp=1,ip
         ipt=(ip-1)*nalt
         jpt=(jp-1)*nalt
         do 8360 ia=1,nalt
         do 8360 ja=1,nalt
            id=ipt+ia
            jd=jpt+ja
            omega(id,jd)=calfa(ia,ja)
            omega(jd,id)=calfa(ia,ja)
8360     continue
8350  continue
c
c     add identity matrix to all diagonal blocks
c
      do 8390 i=1,nt
8390     omega(i,i)=omega(i,i)+1.0
      return
c
c     model=4: calfa on all blocks, plus cnunu matrix on diagonal blocks
c
c     transfer calfa block to all off-diagonal blocks
c
8400  do 8450 ip=1,nper
      do 8450 jp=1,ip-1
         ipt=(ip-1)*nalt
         jpt=(jp-1)*nalt
         do 8460 ia=1,nalt
         do 8460 ja=1,nalt
            id=ipt+ia
            jd=jpt+ja
            omega(id,jd)=calfa(ia,ja)
            omega(jd,id)=calfa(ia,ja)
8460     continue
8450  continue
c
c     add calfa and cnunu matrices on all diagonal blocks
c
      do 8470 ip=1,nper
         ipt=(ip-1)*nalt
         do 8480 ia=1,nalt
         do 8480 ja=1,nalt
            id=ipt+ia
            jd=ipt+ja
            omega(id,jd)=calfa(ia,ja)+cnunu(ia,ja)
8480     continue
8470  continue
      return
c
c     model=7,8: add calfa on each block of omega matrix 
c
8500  do 8510 ia=1,nalt
      do 8510 ja=1,nalt
         caij=calfa(ia,ja)
         do 8520 it=1,nper
         do 8520 jt=1,nper
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
            omega(ita,jta)=omega(ita,jta)+caij
8520     continue
8510  continue
      return
9000  write (*,*) ' over/underflow in omega1'
      if (ispool.ne.0) write (kanal0,*) ' over/underflow in omega1'
      return 1
      end
c
c
      subroutine omega2(parm,np,*)
c=======================================================================
c
c     creates covariance matrix omega and its derivative domega
c     w.r.t. deep parameters
c
c     (c) axel boersch-supan      25. february 1990
c                                 update: 27. august 1990
c
c     multiperiod model:
c       1: iid across periods, iid across alternatives
c       2: iid across periods, correlated across alternatives
c       3: random effects, iid across alternatives
c       4: random effects, correlated across alternatives
c       5: ar1 errors, iid across alternatives
c       6: ar1 errors, correlated across alternatives
c       7: random effects plus ar1 errors, iid across alternatives
c       8: random effects plus ar1 errors, correlated across alt's
c
c=======================================================================
c
c     model:   eps_it = alfa_i + u_it
c
c              where u_it = rho_i*u_it-1 + nu_it
c
c                    var(nu_it)          = sn_i     vec(nalt)
c                    corr(nu_it,nu_jt)   = sn_ij    mat(nalt2)
c
c                    var(alfa_i)         = sa_i     vec(nalt)
c                    corr(alfa_i,alfa_j) = sa_ij    mat(nalt2)
c
c                    corr(u_it,u_it-1)   = rho_i    vec(nalt)
c
c     note reparametrization: rho/corr=(exp(s)-1)/(exp(s)+1)
c                             stddev=exp(s)
c
c
c     offset  end  length            variable
c     ------  ---  ----------------  ----------------------------------
c     1       nx   nx                x-variables
c     nx+1    nb   (ny+nc)*(nalt-1)  y-interactions incl. constants
c     ns1+1   ns2  nalt              std.dev's of unobserved utils
c     ns2+1   ns3  naltud            correlation of unobserved utils
c     ns3+1   ns4  nalt              std.dev's of random effects
c     ns4+1   ns5  naltud            correlation of random effects
c     ns5+1   npp  nalt              autocorrelations
c
c     note: naltud = nalt*(nalt-1)/2
c
c======================================================================
c
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c
      dimension parm(np)
c
      common calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns),
     *       cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns),
     *       omegax(maxnt,maxnt)
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / doutp / kanal0,ispool,ipause
c
c     (a) transfer updated parm vector on full parmp vector
c         and initialize matrices with identity/zero matrix
c
      do 100 i=1,np
100      parmp(loc1(i))=parm(i)
c
      do 110 i=1,nt
         do 111 j=1,nt
            omega(i,j)=0.0
            do 112 k=1,ns
               domega(i,j,k)=0.0
112         continue
111      continue
         omega(i,i)=1.0
110   continue
c
      do 120 i=1,nalt
      do 120 j=1,nalt
      do 120 k=1,ns
         dcalfa(i,j,k)=0.0
         dcnunu(i,j,k)=0.0
120   continue
c
c     (b) correlation matrix of unobserved utility components
c         (cases 2,4)
c
2000  if (model.ne.2.and.model.ne.4) goto 3000
c
      do 210 ia=1,nalt
      do 210 ja=1,nalt
c
c        fetch standard deviation and correlation of unobserved
c        utility components and undo transformation
c
         ns1i=ns1+ia
         ns1j=ns1+ja
         ns2ij=ns2+maps(ia,ja)
         stni=parmp(ns1i)
         if (dabs(stni).gt.707.0) goto 9000
         stni=dexp(stni)
         stnj=parmp(ns1j)
         if (dabs(stnj).gt.707.0) goto 9000
         stnj=dexp(stnj)
         if (ia.eq.ja) then
            cnij=1.0
         else
            temp=parmp(ns2ij)
            if (dabs(temp).gt.707.0) goto 9000
            temp=dexp(temp)
            cnij=(temp-1)/(temp+1)
            dcndp=temp*2/((temp+1)**2)
         end if
c
c        location of derivatives
c
         ms1i =locs(ns1i)
         ms1j =locs(ns1j)
         ms2ij=locs(ns2ij)
c
c        compute covariance of unobserved utilities
c
         cnij=stni*stnj*cnij
         cnunu(ia,ja)=cnij
c
c        derivatives
c
         if (ia.eq.ja) then 
            dcnunu(ia,ja,ms1i )=2*cnij
         else
            dcnunu(ia,ja,ms1i )=cnij
            dcnunu(ia,ja,ms1j )=cnij
            dcnunu(ia,ja,ms2ij)=stni*stnj*dcndp
         end if
210   continue
c
c     (c) correlation matrix of random effects
c         (cases 3,4,7,8)
c
3000  if (model.le.2.or.model.eq.5.or.model.eq.6) goto 4000
c
      do 310 ia=1,nalt
      do 310 ja=1,nalt
c
c        fetch standard deviation and correlation of random effects
c        and undo transformation
c
         ns3i=ns3+ia
         ns3j=ns3+ja
         ns4ij=ns4+maps(ia,ja)
         stai=dexp(parmp(ns3i))
         staj=dexp(parmp(ns3j))
         if (ia.eq.ja) then
            caij=1.0
         else
            temp=dexp(parmp(ns4ij))
            caij=(temp-1)/(temp+1)
            dcadp=temp*2/((temp+1)**2)
         end if
c
c        location of derivatives
c
         ms3i =locs(ns3i)
         ms3j =locs(ns3j)
         ms4ij=locs(ns4ij)
c
c        compute covariance of unobserved utilities
c
         caij=stai*staj*caij
         calfa(ia,ja)=caij
c
c        derivatives
c
         if (ia.eq.ja) then
            dcalfa(ia,ja,ms3i )=2*caij
         else
            dcalfa(ia,ja,ms3i )=caij
            dcalfa(ia,ja,ms3j )=caij
            dcalfa(ia,ja,ms4ij)=stai*staj*dcadp
         end if
310   continue
c
c     (d) autoregressive process with iid nu's
c         (cases 5,7)
c
4000  if (model.ne.5.and.model.ne.7) goto 5000
c
      do 400 ia=1,nalt
      do 400 ja=1,nalt
c
c        fetch rho and undo transformation
c
         ns5i=ns5+ia
         ns5j=ns5+ja
         temp=dexp(parmp(ns5i))
         rhoi=(temp-1)/(temp+1)
         drhidp=2*temp/((temp+1)**2)
         temp=dexp(parmp(ns5j))
         rhoj=(temp-1)/(temp+1)
         drhjdp=2*temp/((temp+1)**2)
c
c        identity for nu's
c
         cnij=0.0
         if (ia.eq.ja) cnij=1.0
c
c        location of derivatives
c
         ms5i =locs(ns5i)
         ms5j =locs(ns5j)
c
c        time dependent part
c
         do 420 it=1,nper
         do 420 jt=1,nper
            is=it-jt
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
            if (jt.lt.it.and.rhoi.ne.0) then
               rhoit=rhoi**is
               ddd=is*rhoit/rhoi
               omega(ita,jta)=rhoit*cnij
c              ---w.r.t. rhoi
               domega(ita,jta,ms5i)=cnij*ddd*drhidp
            else if (it.lt.jt.and.rhoj.ne.0) then
               rhojt=rhoj**(-is)
               ddd=(-is)*rhojt/rhoj
               omega(ita,jta)=rhojt*cnij
c              ---w.r.t. rhoj
               domega(ita,jta,ms5j)=cnij*ddd*drhjdp
            else if (it.eq.jt) then
               omega(ita,jta)=cnij
            end if
            if (rhoi.eq.0.and.is.eq. 1) domega(ita,jta,ms5i)=cnij*drhidp
            if (rhoj.eq.0.and.is.eq.-1) domega(ita,jta,ms5j)=cnij*drhjdp
420      continue
400   continue
c
c     (e) autoregressive process with correlated nu's
c         (cases 6,8)
c
5000  if (model.ne.6.and.model.ne.8) goto 8000
c
      do 500 ia=1,nalt
      do 500 ja=1,nalt
c
c        fetch std.dev's and correlation of unobserved util's
c        and undo transformation
c
         ns1i=ns1+ia
         ns1j=ns1+ja
         ns2ij=ns2+maps(ia,ja)
         stni=dexp(parmp(ns1i))
         stnj=dexp(parmp(ns1j))
         if (ia.eq.ja) then
            cnij=1.0
         else
            temp=dexp(parmp(ns2ij))
            cnij=(temp-1)/(temp+1)
            dcndp=2*temp/((temp+1)**2)
         end if
c
c        fetch rho and and undo transformations
c
         ns5i=ns5+ia
         ns5j=ns5+ja
         temp=dexp(parmp(ns5i))
         rhoi=(temp-1)/(temp+1)
         drhidp=2*temp/((temp+1)**2)
         temp=dexp(parmp(ns5j))
         rhoj=(temp-1)/(temp+1)
         drhjdp=2*temp/((temp+1)**2)
         rho1ij=1.0-rhoi*rhoj
c
c        stationarity assumption
c
         snij=stni*stnj/rho1ij
         cnij=snij*cnij
         if (ia.eq.ja) then
            ddd=2*rhoi/rho1ij
         else
            d1i=rhoi/rho1ij
            d1j=rhoj/rho1ij
         end if
c
c        location of derivatives
c
         ms1i =locs(ns1i)
         ms1j =locs(ns1j)
         ms2ij=locs(ns2ij)
         ms5i =locs(ns5i)
         ms5j =locs(ns5j)
c
c        time dependent part
c
         do 520 it=1,nper
         do 520 jt=1,nper
            is=it-jt
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
c
c           3 cases: it>jt (upper right blocks)
c                    it=jt (diagonal blocks)
c                    it<jt (lower left blocks)
c
            if (jt.lt.it.and.rhoi.ne.0) then
               rhos=rhoi**is
               omega(ita,jta)=rhos*cnij
c              ---w.r.t. std.dev./corr. of unobserved utilities
               if (ia.eq.ja) then
                  domega(ita,jta,ms1i)=2*rhos*cnij
               else
                  domega(ita,jta,ms1i)=rhos*cnij
                  domega(ita,jta,ms1j)=rhos*cnij
                  domega(ita,jta,ms2ij)=rhos*snij*dcndp
               end if
c              ---w.r.t. rhoi,rhoj
               if (ia.eq.ja) then
                  dzz=is/rhoi+ddd
                  domega(ita,jta,ms5i)=rhos*cnij*dzz*drhidp
               else
                  dzz=is/rhoi+d1j
                  domega(ita,jta,ms5i)=rhos*cnij*dzz*drhidp
                  domega(ita,jta,ms5j)=rhos*cnij*d1i*drhjdp
               end if
            else if (it.lt.jt.and.rhoj.ne.0) then
               rhos=rhoj**(-is)
               omega(ita,jta)=rhos*cnij
c              ---w.r.t. std.dev./corr. of unobserved utilities
               if (ia.eq.ja) then
                  domega(ita,jta,ms1i)=2*rhos*cnij
               else
                  domega(ita,jta,ms1i)=rhos*cnij
                  domega(ita,jta,ms1j)=rhos*cnij
                  domega(ita,jta,ms2ij)=rhos*snij*dcndp
               end if
c              ---w.r.t. rhoi,rhoj
               if (ia.eq.ja) then
                  dzz=(-is)/rhoj+ddd
                  domega(ita,jta,ms5i)=rhos*cnij*dzz*drhidp
               else
                  dzz=(-is)/rhoj+d1i
                  domega(ita,jta,ms5i)=rhos*cnij*d1j*drhidp
                  domega(ita,jta,ms5j)=rhos*cnij*dzz*drhjdp
               end if
            else if (it.eq.jt) then
               omega(ita,jta)=cnij
c              ---w.r.t. std.dev./corr. of unobserved utilities
               if (ia.eq.ja) then
                  domega(ita,jta,ms1i)=2*cnij
               else
                  domega(ita,jta,ms1i)=cnij
                  domega(ita,jta,ms1j)=cnij
                  domega(ita,jta,ms2ij)=snij*dcndp
               end if
c              ---w.r.t. rhoi,rhoj
               if (ia.eq.ja) then
                  domega(ita,jta,ms5i)=cnij*ddd*drhidp
               else
                  domega(ita,jta,ms5i)=cnij*d1j*drhidp
                  domega(ita,jta,ms5j)=cnij*d1i*drhjdp
               end if
            end if
            if (rhoi.eq.0.and.is.eq. 1) domega(ita,jta,ms5i)=cnij*drhidp
            if (rhoj.eq.0.and.is.eq.-1) domega(ita,jta,ms5j)=cnij*drhjdp
520      continue
500   continue
c
c     (g) add pieces together, according to case
c
8000  goto (8100,8200,8300,8400,8100,8100,8500,8500),model
c            1    2    3    4    5    6    7    8
c
c     model=1,5,6: done
c
8100  return
c
c     model=2: transfer cnunu block to all diagonal blocks
c
8200  do 8250 ip=1,nper
         ipt=(ip-1)*nalt
         do 8260 ia=1,nalt
         do 8260 ja=1,nalt
            id=ipt+ia
            jd=ipt+ja
            omega(id,jd)=cnunu(ia,ja)
            do 8270 k=1,ns
8270           domega(id,jd,k)=dcnunu(ia,ja,k)
8260     continue
8250  continue
      return
c
c     model=3: calfa on all blocks plus identity matrix on diag. blocks
c
c     transfer calfa block to all blocks
c
8300  do 8350 ip=1,nper
      do 8350 jp=1,ip
         ipt=(ip-1)*nalt
         jpt=(jp-1)*nalt
         do 8360 ia=1,nalt
         do 8360 ja=1,nalt
            id=ipt+ia
            jd=jpt+ja
            omega(id,jd)=calfa(ia,ja)
            omega(jd,id)=calfa(ia,ja)
            do 8370 k=1,ns
               domega(id,jd,k)=dcalfa(ia,ja,k)
8370           domega(jd,id,k)=dcalfa(ia,ja,k)
8360     continue
8350  continue
c
c     add identity matrix to all diagonal blocks
c
      do 8390 i=1,nt
8390     omega(i,i)=omega(i,i)+1.0
      return
c
c     model=4: calfa on all blocks, plus cnunu matrix on diagonal blocks
c
c     transfer calfa block to all off-diagonal blocks
c
8400  do 8450 ip=1,nper
      do 8450 jp=1,ip-1
         ipt=(ip-1)*nalt
         jpt=(jp-1)*nalt
         do 8460 ia=1,nalt
         do 8460 ja=1,nalt
            id=ipt+ia
            jd=jpt+ja
            omega(id,jd)=calfa(ia,ja)
            omega(jd,id)=calfa(ia,ja)
            do 8465 k=1,ns
               domega(id,jd,k)=dcalfa(ia,ja,k)
8465           domega(jd,id,k)=dcalfa(ia,ja,k)
8460     continue
8450  continue
c
c     add calfa and cnunu matrices on all diagonal blocks
c
      do 8470 ip=1,nper
         ipt=(ip-1)*nalt
         do 8480 ia=1,nalt
         do 8480 ja=1,nalt
            id=ipt+ia
            jd=ipt+ja
            omega(id,jd)=calfa(ia,ja)+cnunu(ia,ja)
            do 8485 k=1,ns
8485           domega(id,jd,k)=dcalfa(ia,ja,k)+dcnunu(ia,ja,k)
8480     continue
8470  continue
      return
c
c     model=7,8: add calfa on each block of omega matrix 
c
8500  do 8510 ia=1,nalt
      do 8510 ja=1,nalt
         caij=calfa(ia,ja)
         do 8520 it=1,nper
         do 8520 jt=1,nper
            ita=(it-1)*nalt+ia
            jta=(jt-1)*nalt+ja
            omega(ita,jta)=omega(ita,jta)+caij
            do 8525 k=1,ns
8525           domega(ita,jta,k)=domega(ita,jta,k)+dcalfa(ia,ja,k)
8520     continue
8510  continue
      return
9000  write (*,*) ' over/underflow in omega2'
      if (ispool.ne.0) write (kanal0,*) ' over/underflow in omega2'
      return 1
      end
c
      subroutine clark2(em,dem,sg,dsg,prob,grad,lderiv,*)
c=====================================================================
c
c     clark-approximation: prob (iloop=1) und ableitung grad (iloop=2)
c
c     grad: 1......nb:  deriv's w.r.t. deep parameters in em (=zquer)
c           nb+1...nd:  deriv's w.r.t. deep parameters in sg (=covz)
c
c     copyright (c) axel brsch-supan
c
c     *** adapted to slp: 27. august 1990
c
c=====================================================================
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c
      dimension grad(maxnp),
     *          em(maxnt1),sg(maxnt1,maxnt1),
     *          dem(maxnt1,maxnb),dsg(maxnt1,maxnt1,maxns)
c
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / doutp / kanal0,ispool,ipause
c
      common / dfunc / em1(maxnt1),em2(maxnt1),sg1(maxnt1,maxnt1),
     *                 dem1(maxnt1,maxnp),dem2(maxnt1,maxnp),
     *                 dsg1(maxnt1,maxnt1,maxnp),da(maxnp),dalf(maxnp)
c
c     nb=number of parameters in em=zquer
c     ns=number of parameters in sq=cov(zquer)
c     nt1=dimension of em and zquer =(nalt-1)*nper
c
      nd=nb+ns
c
      if (lderiv) goto 100
c
c     initialize: i=1
c
      do 10 i=1,nt1
         em1(i)=em(i)
         do 10 j=1,nt1
10          sg1(i,j)=sg(i,j)
c
ccccc      write (*,'(a,6f10.4)') 'zquer=',(em(i),i=1,nt1)
cc         write (kanal0,'(a,6f10.4)') 'zquer=',(em(i),i=1,nt1)
cc         mxt1=maxnt1
ccccc      call prtmat('$$$ covz',sg,nt1,nt1,mxt1)
c
c     stepwise approximation: i=2,..,nt-1
c
      do 20 i=2,nt1
         im=i-1
         a      = sg1(im,im) + sg(i,i) - 2.0*sg1(im,i)
         if (a.lt.1.0d-30) goto 99
         a      = dsqrt(a)
         alf    = ( em1(im) - em(i) )/a
         den    = denphi(alf)
         cum    = cumphi(alf)
         em1(i) = em(i) + ( em1(im) - em(i) )*cum + a*den
         em2(i) = em(i)**2 + sg(i,i)
     *           + ( em1(im)**2 + sg1(im,im) - em(i)**2 - sg(i,i) )*cum
     *           + ( em1(im) + em(i) )*a*den
         sg1(i,i) = em2(i) - em1(i)**2
         do 20 j=i+1,nt1
            sg1(i,j) = sg(i,j) + (sg1(im,j) - sg(i,j))*cum
            sg1(j,i) = sg1(i,j)
20    continue
c
c     final probability: i=nt1
c
      a = sg1(nt1,nt1)
      if (a.lt.1.0d-30) goto 99
      a    =  dsqrt(a)
      alf  = -em1(nt1)/a
      prob =  cumphi(alf)
      return
c
c=====version with derivatives:========================
c
c     initialize: i=1
c
100   do 101 i=1,nt1
         em1(i)=em(i)
         do 101 j=1,nt1
101         sg1(i,j)=sg(i,j)
c
      do 104 i=1,nt1
         do 102 k=1,nb
102         dem1(i,k)=dem(i,k)
         do 103 ks=1,ns
            k=nb+ks
            dem (i,k)=0
103         dem1(i,k)=0
104   continue
c
      do 107 i=1,nt1
      do 107 j=1,nt1
         do 105 k=1,nb
105         dsg1(i,j,k)=0
         do 106 ks=1,ns
            k=nb+ks
106         dsg1(i,j,k)=dsg(i,j,ks)
107   continue
c
c     stepwise approximation: i=2,..,nt-1
c
      do 200 i=2,nt1
         im=i-1
c
         a = sg1(im,im) + sg(i,i) - 2.0*sg1(im,i)
         if (a.lt.1.0d-30) goto 99
         a = dsqrt(a)
         do 201 k=1,nb
201         da(k)=0.5*(dsg1(im,im,k)-2.0*dsg1(im,i,k))/a
         do 202 ks=1,ns
            k=nb+ks
            da(k)=0.5*(dsg1(im,im,k)+dsg(i,i,ks)-2.0*dsg1(im,i,k))/a
202      continue
c
         alf = ( em1(im) - em(i) )/a
         do 205 k=1,nd
205         dalf(k)=( dem1(im,k) - dem(i,k) - alf*da(k) )/a
c
         den = denphi(alf)
         cum = cumphi(alf)
c
         em1(i) = em(i) + ( em1(im) - em(i) )*cum + a*den
         do 210 k=1,nd
            dem1(i,k) = dem(i,k)
     *                + ( dem1(im,k) - dem(i,k) )*cum
     *                + ( em1(im) - em(i) )*den*dalf(k)
     *                + da(k)*den + a*den*(-alf)*dalf(k)
210      continue
c
         em2(i) = em(i)**2 + sg(i,i)
     *           + ( em1(im)**2 + sg1(im,im) - em(i)**2 - sg(i,i) )*cum
     *           + ( em1(im) + em(i) )*a*den
         do 215 k=1,nb
            dem2(i,k) = 2.0*em(i)*dem(i,k)
     *      + (  2.0*em1(im)*dem1(im,k) + dsg1(im,im,k)
     *          -2.0*em (i )*dem (i ,k)                  )*cum
     *      + ( em1(im)**2+sg1(im,im)-em(i)**2-sg(i,i) )*den*dalf(k)
     *      + (dem1(im,k)+dem(i,k))*a*den
     *      + (em1(im)+em(i))*da(k)*den
     *      + (em1(im)+em(i))*a*den*(-alf)*dalf(k)
215      continue
         do 216 ks=1,ns
            k=nb+ks
            dem2(i,k) = 2.0*em(i)*dem(i,k) + dsg(i,i,ks)
     *      + (  2.0*em1(im)*dem1(im,k) + dsg1(im,im,k )
     *          -2.0*em (i )*dem (i ,k) - dsg (i ,i ,ks) )*cum
     *      + ( em1(im)**2+sg1(im,im)-em(i)**2-sg(i,i) )*den*dalf(k)
     *      + (dem1(im,k)+dem(i,k))*a*den
     *      + (em1(im)+em(i))*da(k)*den
     *      + (em1(im)+em(i))*a*den*(-alf)*dalf(k)
216      continue
c
         sg1(i,i) = em2(i) - em1(i)**2
         do 220 k=1,nd
            dsg1(i,i,k) = dem2(i,k) - 2.0*em1(i)*dem1(i,k)
220      continue
c
         do 230 j=i+1,nt1
            sg1(i,j) = sg(i,j) + (sg1(im,j) - sg(i,j))*cum
            sg1(j,i) = sg1(i,j)
            do 231 k=1,nb
               dsg1(i,j,k) = ( dsg1(im,j,k)              )*cum
     *                     + ( sg1(im,j) - sg(i,j) )*den*dalf(k)
               dsg1(j,i,k) = dsg1(i,j,k)
231         continue
            do 232 ks=1,ns
               k=nb+ks
               dsg1(i,j,k) = dsg(i,j,ks)
     *                     + ( dsg1(im,j,k) - dsg(i,j,ks) )*cum
     *                     + ( sg1(im,j) - sg(i,j) )*den*dalf(k)
               dsg1(j,i,k) = dsg1(i,j,k)
232         continue
230      continue
c
200   continue
c
c     final probability: i=nt1
c
      a = sg1(nt1,nt1)
      if (a.lt.1.0d-30) goto 99
      a = dsqrt(a)
      do 301 k=1,nd
301      da(k)=0.5*dsg1(nt1,nt1,k)/a
      alf = -em1(nt1)/a
      do 302 k=1,nd
302      dalf(k)=-( dem1(nt1,k) + alf*da(k) )/a
      prob =  cumphi(alf)
      den  =  denphi(alf)
      do 303 k=1,nd
303      grad(k)=den*dalf(k)
c
      return
99    write (*,'(a,i2,d15.7)')
     *       ' clarkk: sqrts0 in step,loop ',nt1,a
      if (ispool.ne.0) write (kanal0,'(a,i2,d15.7)')
     *       ' clarkk: sqrts0 in step,loop ',nt1,a
      return 1
      end
c
c
      subroutine logit2(zquer,dzquer,prob,grad,lderiv,*)
c=====================================================================
c
c     multinomial logit model (after ialt transformation)
c
c     *** adapted to sml: 27. august 1990
c
c=====================================================================
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
      dimension zquer(maxnt1),dzquer(maxnt1,maxnb),grad(maxnp)
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dfunc / ezq(maxnt1)
c
      prob=1.0
      if (lderiv) then
         do 5 k=1,nb
5           grad(k)=0.0
      end if
c
      mp0=0
      do 100 iper=1,nper
c
         dinc=1.d0
         do 10 i=1,nalt1
            it1=mp0+i
            zq=zquer(it1)
            if (zq.gt.500.d0.or.zq.lt.-500.d0) then
cover          write (*,*) 'under/overflow in logit: zq=',zq
               return 1
            end if
            ezq(i)=dexp(zq)
            dinc=dinc+ezq(i)
10       continue
c
         if (dinc.gt.1.d125.or.dinc.lt.-1.d125) then
cover       write (*,*) 'under/overflow in logit: dinc=',dinc
            return 1
         end if
c
         prob1=1.d0/dinc
         prob0=prob
         prob=prob*prob1
c
         if (lderiv) then
            do 20 k=1,nb
               sum=0.d0
               do 25 j=1,nalt1
                  jt1=mp0+j
25                sum=sum+ezq(j)*dzquer(jt1,k)
               grad1=-sum/dinc**2
               grad(k)=grad(k)*prob1+prob0*grad1
20          continue
         end if
c
      mp0=mp0+nalt1
100   continue
c
      return
      end
c-----------------------------------------------------------------------
c.pl100/.rm72/
c     sml-inp: input routines for mpp-msl programm
c
c     (c) axel boersch-supan,  1. september 1990
c
c----------------------------------------------------------------------
      subroutine inparm(parm,np,*)
c-----------------------------------------------------------------------
c
c     profile and parameter input for mpp-msl
c
c     two parameter vectors are involved:
c       parmp: vector of primary x and y parameters (length=nb+ns)
c              this vector has a fixed structure with zeroes for
c              parameters not involved in estimation
c       parm:  complete vector of estimated parameters (length=np)
c
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
c
      dimension   parm(maxnp)
      character*8 pname(maxnpp)
      character*20 runtag
      integer*4   iseed
      integer*2   iarg1,ierr1
      character*40 text,filen,spfile,dafile,pafile
cvah:start
      character*80 cline
      integer*4 istrt(5),iend(5),narg
cvah:end
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / dname / pname
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / doutp / kanal0,ispool,ipause
      common / doptm / mchk0,mchk1,iterl1,acc1,ialg1,
     *                       mchk2,iterl2,acc2,ialg2,nwesml
      common / dsimu / nsimul
      common / dseed / iseed
      common / dtag  / runtag
      common / dfile / spfile,pafile
c
c     initialize before reading profile
c
      kanal0=10
      ispool=0
      ipause=-1
c
      kanal=20
      text='profile/parameter-input-file'
      mode=1
      iarg1=1
caxel      call getarg(iarg1,filen,ierr1)
cvah:start
      call parsecl(cline,istrt,iend,narg)
      ierr1=0
      if(narg.gt.0) then
         filen=cline(istrt(1):iend(1))
         ierr1=1
      endif
cvah:end
      if (ierr1.le.0) then
         write (*,*) 'enter ',text,'[def=sml.pro]'
         read (*,'(a40)') filen
         if (filen.eq.' ') filen='sml.pro'
      end if
      call openf(kanal,text,filen,mode)
c
c======================================================================
c
c *** read in profile and open files ******************************
c
c     ipause = seconds pause between screens
c
c     spfile = file name for spool output
c     ispool = indicator for spool output
c
c     dafile = file name for data input
c     pafile = file name for parameter output
c
c     mchk0-2= derivative checks for omega, initial est., sml estimation
c     iterl  = number of iterations of optimizer
c     acc 1-2= accuracy of optimizer                 for init, sml
c     ialg   = optimizer (2=dfp, 4=bhhh)
c
c     method = method for initial estimator (0=none, 1=logit, 2=clark)
c     nsimul = number of replications in simulation
c     iseed  = first seed for simulation
c
c======================================================================
c
      m1=1
      m2=2
      m3=3
      call skip(kanal,m2)
      read (kanal,*) ipause
c
      read (kanal,'(/3x,a40)') spfile
      ispool=0
      if (spfile.ne.' ') then
         text='spool-file'
         kanal0=10
         mode=3
         call openf(kanal0,text,spfile,mode)
         ispool=1
         write (kanal0,1) runtag,
     *                    maxalt,maxper,maxnt,maxnt1,
     *                    maxnb,maxns,maxnp,maxnpp,
     *                    maxbuf
1        format(
     *//' ============================================================'
     * /'       **** multinomial multiperiod probit model ****'
     * /'    version',        a20     ,'  (c) axel boersch-supan'
     * /' ============================================================'
     * /'       maxalt=',i3,' maxper=',i3,' maxnt=',i3,' maxnt1=',i3
     * /'       maxnb =',i3,' maxns =',i3,' maxnp=',i3,' maxnpp=',i3
     * /'                                         maxbuf=',i6
     * /' ============================================================')
      end if
      read (kanal,'(/3x,a40)') dafile
        text='data-file'
        kanal7=7
        mode=1
        call openf(kanal7,text,dafile,mode)
      read (kanal,'(/3x,a40)') pafile
        runtag='                    '
        runtag=pafile(1:20)
        text='result-file'
        kanal8=8
        mode=3
        call openf(kanal8,text,pafile,mode)
        write (kanal8,*) 'sml output file ',runtag
        close(kanal8)
c
      call skip(kanal,m2)
      read (kanal,*) mchk0
      call skip(kanal,m2)
      read (kanal,*) method,      mchk1,iterl1,acc1,ialg1
      call skip(kanal,m3)
      read (kanal,*) nsimul,iseed,mchk2,iterl2,acc2,ialg2
      call skip(kanal,m2)
      read (kanal,*) nwesml
c
      write (*,10)
     *      ispool,ipause,method,nsimul,iseed,mchk0,mchk1,mchk2,nwesml
      if (ispool.ne.0) write (kanal0,10)
     *      ispool,ipause,method,nsimul,iseed,mchk0,mchk1,mchk2,nwesml
10    format(///' multinomial-multiperiod probit model'
     *         /' ===================================='
     *         /' current profile:'
     *         /' -----------------------------'
     *         /' spool file activated:  ', i6
     *         /' pause/batch mode:      ', i6
     *         /' initial estimator:     ', i6
     *         /' number of replications:', i6
     *         /' initial seed:      ', i10
     *         /' derivative checks:     '
     *         /' - omega                ', i6
     *         /' - initial estimator    ', i6
     *         /' - sml estimator        ', i6
     *         /' robust covariance:     ', i6)
      call snooze(ipause)
c
c======================================================================
c
c *** read in dimensions *************************************
c
c     nindiv= number of individuals (default)
c     nalt  = number of alternatives, indexed i=1,...,nalt in y
c     nper  = largest number of waves in panel
c     nx    = number of alternative-specific explanatory variables
c     ny    = number of agent-specific explanatory variables
c     nc    = alternative-specific constants present
c     nwt   = weights present
c
c======================================================================
c
      call skip(kanal,m2)
      read (kanal,*) nindiv
      call skip(kanal,m1)
      read (kanal,*) nalt,nx,ny,nc,nwt
      nalt1=nalt-1
      nalt2=nalt-2
      naltc=nalt1*nalt2/2
      naltud=nalt*nalt1/2
      if (nc.gt.0) nc=1
      if (nc.lt.0) nc=0
      if (nwt.gt.0) nwt=1
      if (nwt.lt.0) nwt=0
      nxy=nx+nalt1*(ny+nc)
      ndd=nx*nalt
      ndw=nx*nalt+1+nwt
      nda=nx*nalt+1+nwt+ny
c
      write (*,11) nalt,nx,ny,nc,nwt,nxy
      if (ispool.ne.0) write (kanal0,11) nalt,nx,ny,nc,nwt,nxy
11    format(///' multinomial-multiperiod probit model'
     *         /' ===================================='
     *         /' number of alternatives:', i6
     *         /' number of x-variables: ', i6
     *         /' number of y-variables: ', i6
     *         /' constants included:    ', i6
     *         /' weights included:      ', i6
     *        //' total number of betas: ', i6)
      call snooze(ipause)
      if (nalt.gt.maxalt) goto 9003
c 
c *** read in error process ********************************************                         
c                                                            ran ar1 mnp
c     multiperiod model:
c     1: iid across periods, iid across alternatives           0   0   0
c     2: iid across periods, correlated across alternatives    0   0   1                                                                     
c     3: random effects, iid across alternatives               1   0   0
c     4: random effects, correlated across alternatives        1   0   1
c     5: ar1 errors, iid across alternatives                   0   1   0
c     6: ar1 errors, correlated across alternatives            0   1   1
c     7: random effects + ar1 errors, iid across alt's         1   1   0
c     8: random effects + ar1 errors, correlated across alt's  1   1   1
c
c=======================================================================
c
      call skip(kanal,m1)
      read (kanal,*) nper,modran,modar1,modmnp
20    if (method.eq.1) then
         model=1
         modran=0
         modar1=0
         modmnp=0
      end if
      if (nper.eq.1) then
         model=1
         if (modmnp.eq.1) model=2
         modran=0
         modar1=0
      else
         model=1
         if (modmnp.eq.1) model=2
         if (modran.eq.1) model=model+2
         if (modar1.eq.1) model=model+4
      end if
c
      ns1=nxy
      ns2=ns1+nalt
      ns3=ns2+naltud
      ns4=ns3+nalt
      ns5=ns4+naltud
      npp=ns5+nalt
c
      nt=nper*nalt
      nt1=nper*nalt1
      ndp=nda*nper
c
      write (*,22) model,nper,modmnp,modran,modar1,npp
      if (ispool.ne.0) write (kanal0,22)
     *             model,nper,modmnp,modran,modar1,npp
22    format(///' specification of error process (model',i2,')'
     *         /' ========================================'
     *         /' number of periods:       ',i2
     *         /' inter-alternative corr:   ',i1
     *         /' random effects:           ',i1
     *         /' autocorrelation:          ',i1
     *         //' required initial values: ',i2)
      call snooze(ipause)
      if (nper.gt.maxper) goto 9005
      if (nt  .gt.maxnt ) goto 9006
      if (npp .gt.maxnpp) goto 9007
c
      if (ipause.ge.0) then
         write (*,*) 'hit enter to accept model, 1 if not>'
         read (*,'(i1)') iquery
         if (iquery.ne.0) then
            write (*,*) 'enter nper, mod_mnp, mod_ran, mod_ar1 >'
            read (*,*) nper,modmnp,modran,modar1
            goto 20
         end if
      end if
c
28    nbuff=nindiv*ndp
      if (nbuff.gt.maxbuf) then
         write (*,29) nindiv,ndp,maxbuf
         if (ispool.ne.0) write (kanal0,29) nindiv,ndp,maxbuf
29       format(' *** error: out of memory ***'
     *         /' nindiv=',i9,' at ndp=',i4,' exceeds maxbuf=',i9)
         if (ipause.lt.0) return 1
         write (*,*) 'enter new nindiv'
         read (*,*) nindiv
         goto 28
      end if
c
c======================================================================
c
c *** read in initialvalues ****************************************
c
c     parmp = primary parameter vector with length npp
c     parm  = estimated parameter vector with length np
c
c     pname = label of variable
c
c     loc   = 0 if fixed
c             1 if estimated
c
c     note: the specification of model overrides istat.
c
c     npp = number of primary parameters (any istat)
c     np  = number of estimated parameters (istat=0)
c
c
c----------------------------------------------------------------------
c
c     number of....
c     estimated   all parametrs:           contents:
c     ---------   --------------------     --------
c     nb          nxy=nx+nalt1*(ny+nc)     beta in xbeta
c     ns          3*nalt+2*naltud          sigma in omega
c     np          npp                      stacked(beta,sigma)
c     ---------   --------------------     --------
c     parm        parmp
c
c
c----------------------------------------------------------------------
c
c     loc    first status (1=estim, 0=fixed), then
c            location of i-th full parameter in estimated parm,
c                    if parameter is estimated, 0 if fixed
c            e.g.:
c            grad(loc(ip)) is derivative associated with
c                          full parameter component parm(ip)
c
c     loc1   location of i-th estimated parameter in full parmp,
c            e.g.:
c            pname(loc1(ip1)) is parameter name associated with
c                          estimated parameter parm(ip1)
c
c----------------------------------------------------------------------
c
c     parmp:
c     offset  end  length            variable
c     ------  ---  ----------------  ----------------------------------
c     1       nx   nx                x-variables
c     nx+1    nxy  (ny+nc)*nalt      y-interactions incl. constants
c     ns1+1   ns2  nalt              std.dev's of unobserved utils
c     ns2+1   ns3  naltud            correlation of unobserved utils
c     ns3+1   ns4  nalt              std.dev's of random effects
c     ns4+1   ns5  naltud            correlation of random effects
c     ns5+1   npp  nalt              autocorrelations
c
c     notes: naltud = nalt*(nalt-1)/2
c
c======================================================================
c
300   m13=13
      call skip(kanal,m13)
      write (*,301)
      if (ispool.ne.0) write (kanal0,301)
301   format(//' initialvalues of parameters:'
     *        /' ============================'/)
      do 310 k=1,npp
         read (kanal,311,end=997) pname(k),loc(k),parmp(k)
311      format(a8,i2,f15.9)
         write (*,312) pname(k),loc(k),parmp(k)
         if (ispool.ne.0)
     *      write (kanal0,312) pname(k),loc(k),parmp(k)
312      format(1x,a8,1x,i2,1x,f15.9)
         if (k.eq.15.or.k.eq.30) call snooze(ipause)
310   continue
      if (npp.ne.15.and.npp.ne.30) call snooze(ipause)
      close(kanal)
c
c     impose model structure.  obey the following restrictions:
c
c     nalt=4:   | var1 cov1 cov2  0 |    - restrictions for i,j=nalt
c               |      var2 cov3  0 |           due to differencing
c               |              1  0 |    - restriction for var3 due to
c               |    0    0    0  1 |           identification of beta
c
      np=0
      nb=0
      ns=0
c
c     --- x and y var's
      do 320 i=1,nxy
         if (loc(i).ne.0) then
            np=np+1
            nb=nb+1
            loc1(np)=i
            loc(i)=np
            end if
320      continue
c
c     ---interalternative correlations
      if (modmnp.eq.1) then
         do 330 i=1,nalt
            if (loc(ns1+i).ne.0) then
               np=np+1
               ns=ns+1
               loc1(np)=ns1+i
               loc(ns1+i)=np
            end if
330      continue
         do 331 i=1,naltud
            if (loc(ns2+i).ne.0) then
               np=np+1
               ns=ns+1
               loc1(np)=ns2+i
               loc(ns2+i)=np
            end if
331      continue
      else
         do 335 i=1,nalt
335         loc(ns1+i)=0
         do 336 i=1,naltud
336         loc(ns2+i)=0
      end if
c
c     ---random effects
      if (modran.eq.1) then
         do 340 i=1,nalt
            if (loc(ns3+i).ne.0) then
               np=np+1
               ns=ns+1
               loc1(np)=ns3+i
               loc(ns3+i)=np
            end if
340      continue
         do 341 i=1,naltud
            if (loc(ns4+i).ne.0) then
               np=np+1
               ns=ns+1
               loc1(np)=ns4+i
               loc(ns4+i)=np
            end if
341      continue
      else
         do 345 i=1,nalt
345         loc(ns3+i)=0
         do 346 i=1,naltud
346         loc(ns4+i)=0
      end if
c
      if (modar1.eq.1) then
         do 350 i=1,nalt
            if (loc(ns5+i).ne.0) then
               np=np+1
               ns=ns+1
               loc1(np)=ns5+i
               loc(ns5+i)=np
            end if
350      continue
      else
         do 355 i=1,nalt
355         loc(ns5+i)=0
      end if
c
      write (*,360) npp,nxy,npp-nxy,np,nb,ns
      if (ispool.ne.0) write (kanal0,360) npp,nxy,npp-nxy,np,nb,ns
360   format(//' number of parameters:'
     *        /' =========================='
     *        /' primary parameters:  ',i3
     *        /'    of those: beta      ',i3
     *        /'    of those: sigma     ',i3
     *        /' estimated parameters:',i3
     *        /'    of those: beta      ',i3
     *        /'    of those: sigma     ',i3)
      call snooze(ipause)
      if (np.eq.0) goto 996
      if (nb.gt.maxnb) goto 9008
      if (ns.gt.maxns) goto 9009
c
c     maps: location of correlation-matrix elements in parmp
c
      k=0
      do 380 i=1,nalt
         maps(i,i)=0
         do 385 j=i+1,nalt
            k=k+1
            maps(i,j)=k
            maps(j,i)=k
385      continue
380   continue
c     write (*,*) 'maps='
c     do 389 i=1,nalt
c        write (*,'(10i3)') (maps(i,j),j=1,nalt)
389   continue
c     call snooze(ipause)
c
c     locs: storage of derivatives in domega
c
      do 390 i=1,npp
390      locs(i)=max(0,loc(i)-nb)
c
c     reparametrize and bound away from -1,0,1, respectively
c
      eps= 0.0001d0
      elim=0.9999d0
      do 401 i=1,nalt
c        ---- std.dev nu
         ppp=parmp(ns1+i)
         if (ppp.lt.0) goto 996
         if (ppp.lt.eps) ppp=eps
         parmp(ns1+i)=log(ppp)
c        ---- std.dev alfa
         ppp=parmp(ns3+i)
         if (ppp.lt.0) goto 996
         if (ppp.lt.eps) ppp=eps
         parmp(ns3+i)=log(ppp)
c        ---- rho
         ppp=parmp(ns5+i)
         if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996
         if (ppp.lt.-elim) ppp=-elim
         if (ppp.gt. elim) ppp= elim
         parmp(ns5+i)=log((1+ppp)/(1-ppp))
401   continue
      do 402 i=1,naltud
c        ---- correl nu
         ppp=parmp(ns2+i)
         if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996
         if (ppp.lt.-elim) ppp=-elim
         if (ppp.gt. elim) ppp= elim
         parmp(ns2+i)=log((1+ppp)/(1-ppp))
c        ---- correl alfa
         ppp=parmp(ns4+i)
         if (ppp.lt.-1.0.or.ppp.gt.1.0) goto 996
         if (ppp.lt.-elim) ppp=-elim
         if (ppp.gt. elim) ppp= elim
         parmp(ns4+i)=log((1+ppp)/(1-ppp))
402   continue
c
c     finally: echo and transfer on actual parameter vector
c
500   write (*,501)
      if (ispool.ne.0) write (kanal0,501)
501   format(//' parameters after transformation:'
     *        /' ================================'/)
      do 510 k=1,npp
         write (*,312) pname(k),loc(k),parmp(k)
         if (ispool.ne.0)
     *      write (kanal0,312) pname(k),loc(k),parmp(k)
         if (k.eq.15.or.k.eq.30) call snooze(ipause)
510   continue
      if (npp.ne.15.and.npp.ne.30) call snooze(ipause)
      do 512 i=1,np
         parm(i)=parmp(loc1(i))
512   continue
c
c     exit
c
900   close (kanal)
      return
996   write (*,*) '*** error: initial values not admissible ***'
      if (ispool.ne.0) write (kanal0,*)
     *            '*** error: initial values not admissible ***'
      return 1
997   write (*,*) '*** error: not enough initial values provided ***'
      if (ispool.ne.0) write (kanal0,*)
     *            '*** error: not enough initial values provided ***'
      return 1
9003  write (*,9013) nalt,maxalt
      if (ispool.ne.0) write (kanal0,9013) nalt,maxalt
9013  format(' *** error: out of memory ***'
     *      /' nalt=',i9,' exceeds max=',i9)
      return 1
9005  write (*,9015) nper,maxper
      if (ispool.ne.0) write (kanal0,9015) nper,maxper
9015  format(' *** error: out of memory ***'
     *      /' nper=',i9,' exceeds max=',i9)
      return 1
9006  write (*,9016) nt,maxnt
      if (ispool.ne.0) write (kanal0,9016) nt,maxnt
9016  format(' *** error: out of memory ***'
     *      /' nt=nper*nalt=',i9,' exceeds max=',i9)
      return 1
9007  write (*,9017) npp,maxnp
      if (ispool.ne.0) write (kanal0,9017) npp,maxnp
9017  format(' *** error: out of memory ***'
     *      /' npp=',i9,' exceeds max=',i9)
      return 1
9008  write (*,9018) nb,maxnb
      if (ispool.ne.0) write (kanal0,9018) nb,maxnb
9018  format(' *** error: out of memory ***'
     *      /' nb=',i9,' exceeds max=',i9)
      return 1
9009  write (*,9019) ns,maxns
      if (ispool.ne.0) write (kanal0,9019) ns,maxns
9019  format(' *** error: out of memory ***'
     *      /' ns=',i9,' exceeds max=',i9)
      return 1
      end
c
      subroutine skip(kanal,lines)
      implicit integer (i-n)
      do 10 i=1,lines
10       read (kanal,'(a1)')
      return
      end
c
      subroutine indata(*)
c-----------------------------------------------------------------------
c
c     initial data input for mpp-msl
c     - reads formatted "indata" file
c     - checks data for outliers (>10.000?)
c     - calculates weighted and unweighted means
c
c-----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
      real*4 data(maxbuf)
      character*72 fmtdat
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / ddata / data
      common / dfreq / wtsum,share(maxnt),wshare(maxnt),
     *                       amean(maxndm),wmean(maxndm)
      common / doutp / kanal0,ispool,ipause
c
c *** dialog (note: file have been opended in inparm) ***
c
      kanal7=7
      read (kanal7,'(a72)') fmtdat
c
      write (*,10) nindiv,nper
      if (ispool.ne.0) write (kanal0,10) nindiv,nper
10    format(//' now starting to read ',i6,' individuals at nper=',i3
     *        /' ===================================================')
      nread=0
      if (ipause.ge.0) then
         write (*,11)
11       format(' enter new number of individuals to read  [0=keep] >')
         read (*,'(bn,i4)') nread
      end if
      if (nread.gt.nindiv) then
         write (*,*) 
     *       'note: truncated to maximum memory allocation'
         if (ispool.ne.0) write (kanal0,*)
     *       'note: truncated to maximum memory allocation'
         nread=0
         call snooze(ipause)
      end if
      if (nread.le.0) nread=nindiv
c
c *** offsets and lengths
c
      ndy=1+nwt+ny
c     ---just to remember:
c     ndd=nx*nalt
c     nda=nx*nalt+1+nwt+ny
c     ---storage for means:
      ndm=ndd+ny
c
c *** initializations for means and frequencies
c
      do 20 i=1,ndm*nper
         amean(i)=0.0
20       wmean(i)=0.0
      do 21 i=1,nalt*nper
         share(i)=0.0
21       wshare(i)=0.0
c
c *** loop through data
c
      nn=0
      iwt=0
      nobs=0
      nindiv=0
      do 100 ihhd=1,nread
      do 101 iper=1,nper
c        ---check eof before all periods read
         ipe1=iper
         if (nn.gt.maxbuf) goto 993
c        ---read data according to whether alt-spec var's exist
         if (nx.gt.0) then
            read (kanal7,fmtdat,err=298,end=299)
     *                    (data(nn+(j-1)*nalt+1),j=1,nx),
     *                    (data(nn+ndd+j),j=1,ndy)
            do 110 i=2,nalt
               read (kanal7,fmtdat,err=991,end=991)
     *                    (data(nn+(j-1)*nalt+i),j=1,nx)
110         continue
         else
            read (kanal7,fmtdat,err=298,end=299)
     *                    (data(nn+ndd+j),j=1,ndy)
         end if
c        ---counter
         nobs=nobs+1
c        ---dependent variable in choice set?
         ic=data(nn+ndd+1)
         if (ic.lt.1.or.ic.gt.nalt) then
            write (*,121) nobs,ic
            if (ispool.ne.0) write (kanal0,121) nobs,ic
121         format(' invalid chosen alternative in observation',i6,i3)
            goto 999
         end if
c        ---extract relevant weight
         wt=1.0
         if (nwt.eq.1) wt=data(nn+ndd+2)
         if (wt.ne.1.0) iwt=1
         if (wt.lt.0) then
            write (*,122) nobs,wt
            if (ispool.ne.0) write (kanal0,122) nobs,wt
122         format(' invalid weight in observation',i6,g20.10)
            goto 999
         end if
c        ---accumulate sample frequencies
         icp=(iper-1)*nalt+ic
         share(icp)  = share(icp)  + 1.0
         wshare(icp) = wshare(icp) + wt
c        ---check alt-specific attributes for typos and accumulate means
         indn=0
         indm=(iper-1)*ndm
         do 140 k=1,nx
         do 140 i=1,nalt
            indn=indn+1
            indm=indm+1
            da=data(nn+indn)
            if (abs(da).gt.10000.0) then
               write (*,142) i,k,nobs,da
               if (ispool.ne.0) write (kanal0,142) i,k,nobs,da
142            format(' data check: x(ialt,k,iobs)=',3i5,g20.10)
            end if
            amean(indm) = amean(indm) + da
            wmean(indm) = wmean(indm) + da*wt
140      continue
c        ---check agent-specific characteristics and accumulate means
         indn=ndd+1+nwt
         do 145 k=1,ny
            indn=indn+1
            indm=indm+1
            da=data(nn+indn)
            if (abs(da).gt.10000.0) then
               write (*,146) k,nobs,da
               if (ispool.ne.0) write (kanal0,146) k,nobs,da
146            format(' data check: y(k,iobs)=',2i5,g20.10)
            end if
            amean(indm) = amean(indm) + da
            wmean(indm) = wmean(indm) + da*wt
145      continue
c        ---new offset
         nn=nn+nda
101   continue
      nindiv=nindiv+1
100   continue
      goto 300
c
c *** end of file
c
298   write (*,*) 'warning: err condition ended reading data'
      write (*,*) '         check format for data file'
299   if (nobs.eq.0) then
         write (*,*) ' *** error: data-file empty ***'
         return 1
      end if
      close(kanal7)
      if (ipe1.ne.1) goto 991
c
c *** sample means and frequencies
c
300   wtsum=0.d0
      do 320 i=1,nalt*nper
320      wtsum=wtsum+wshare(i)
      do 321 i=1,ndm*nper
         amean(i)=amean(i)/nindiv
321      wmean(i)=wmean(i)/wtsum*nper
      do 322 i=1,nalt*nper
         share(i)=share(i)/nindiv
322      wshare(i)=wshare(i)/wtsum*nper
c
400   write (*,404) nindiv,nobs,wtsum
      if (ispool.ne.0) write (kanal0,404) nindiv,nobs,wtsum
404   format(//' unweighted number of individuals: ',i7
     *       /' unweighted number of observations:',i7
     *       /' weighted number of observations:',f9.1
     *       /' =========================================')
c
      write (*,409) 'unweighted sample frequencies by period:'
      write (*,411) (i,i=1,nalt)
      do 410 iper=1,nper
410      write (*,412) iper,(share((iper-1)*nalt+i),i=1,nalt)
      call snooze(ipause)
      write (*,409) 'unweighted sample averages by period:'
      write (*,421) (i,i=1,ndm)
      do 420 iper=1,nper
420      write (*,422) iper,(amean((iper-1)*ndm+i),i=1,ndm)
      call snooze(ipause)
c
409   format(/' ',a/' -----------------------------------------')
411   format(' ialt='    ,10i7)
412   format(' t=',i2,':',10f7.3)
421   format(' ivar='    ,10i7:  ,5(/6x,10i7:))
422   format(' t=',i2,':',10f7.3:,5(/6x,10f7.3:))
c
      if (iwt.eq.1) then
         write (*,409) 'weighted sample frequencies by period:'
         write (*,411) (i,i=1,nalt)
         do 415 iper=1,nper
415         write (*,412) iper,(wshare((iper-1)*nalt+i),i=1,nalt)
         call snooze(ipause)
         write (*,409) 'weighted sample averages by period:'
         write (*,421) (i,i=1,ndm)
         do 425 iper=1,nper
425         write (*,422) iper,(wmean((iper-1)*ndm+i),i=1,ndm)
         call snooze(ipause)
      else
         write (*,409) 'note: all weights are one'
         call snooze(ipause)
      end if
c
      if (ispool.eq.0) goto 900
      k=kanal0
      write (k,409) 'unweighted sample frequencies by period:'
      write (k,411) (i,i=1,nalt)
      do 414 iper=1,nper
414      write (k,412) iper,(share((iper-1)*nalt+i),i=1,nalt)
      write (k,409) 'unweighted sample averages by period:'
      write (k,421) (i,i=1,ndm)
      do 424 iper=1,nper
424      write (k,422) iper,(amean((iper-1)*ndm+i),i=1,ndm)
c
      if (iwt.eq.1) then
         write (k,409) 'weighted sample frequencies by period:'
         write (k,411) (i,i=1,nalt)
         do 419 iper=1,nper
419         write (k,412) iper,(wshare((iper-1)*nalt+i),i=1,nalt)
         write (k,409) 'weighted sample averages by period:'
         write (k,421) (i,i=1,ndm)
         do 429 iper=1,nper
429         write (k,422) iper,(wmean((iper-1)*ndm+i),i=1,ndm)
      else
         write (k,409) 'note: all weights are one'
      end if
c
900   return
991   write (*,992)
992   format(/' alignment error:'/' at least one observation',
     *        ' does not have records for all alternatives.')
999   write (*,*) '*** error in input data set ***'
      return 1
993   write (*,994)
994   format(/' *** fatal error: size of data set exceeds memory ***')
      return 1
      end
c-----------------------------------------------------------------------
c.pl100/.rm72/
c     sml-util: external utilities for method of simulated likelihood
c
c     copyright (c) axel boersch-supan
c     version  26. february 1990, update from mpputil 1. july 1990
c                                 update from mslutil 27. august 1990
c
c     subroutines: dcomeg  derivative check of omega
c                  fdiff   numerical first derivatives
c                  sdiff   numerical second derivatives
c                  promeg  print omega
c                  prtmat  print matrix
c                  mslout  result file
c                  cwesml  wesml-robust covariance computation
c
c----------------------------------------------------------------------
      subroutine dcomeg(parm,np,eps)
c----------------------------------------------------------------------
c
c     checks derivatives of omega with respect to deep sigma parameters
c
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
c
      dimension parm(np)
      character*40 spfile,pafile
c
c     common / domeg1/
      common
     *       calfa(maxalt,maxalt),dcalfa(maxalt,maxalt,0:maxns),
     *       cnunu(maxalt,maxalt),dcnunu(maxalt,maxalt,0:maxns),
     *       omegax(maxnt,maxnt)
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / doutp / kanal0,ispool,ipause
      common / dfile / spfile,pafile
c
c     analytic derivs
c
      write (*,*) 'derivative check at eps=',eps
      call omega2(parm,np,*99)
      write (*,'(/a)') ' omega matrix='
      if (ispool.ne.0) write (kanal0,'(/a)') ' omega matrix='
      do 215 i=1,nt
         write (*,'(12f6.3))') (omega(i,j),j=1,nt)
         if (ispool.ne.0) write (kanal0,'(12f6.3))') (omega(i,j),j=1,nt)
         do 215 j=1,nt
215         omegax(i,j)=omega(i,j)
      if (ipause.ge.0) pause 'hit enter for more....'
c
c     numerical derivs
c
      do 220 k=1,ns
         parm(nb+k)=parm(nb+k)+eps
         call omega1(parm,np,*99)
         devmax=-1.0
         do 221 i=1,nt
         do 221 j=1,nt
            omega(i,j)=(omega(i,j)-omegax(i,j))/eps
            devij=dabs(omega(i,j)-domega(i,j,k))
            if (devij.gt.devmax) devmax=devij
221      continue
         write (*,222) k,devmax
         if (ispool.ne.0) write (kanal0,222) k,devmax
222      format(//' analytical/numerical derivs w.r.t. s',i2.2,
     *            '  ***  max abs error=',d10.3,'  ***'/' ',74('='))
         do 224 i=1,nt
            write (*,223) ' a:',(domega(i,j,k),j=1,nt)
            write (*,223) ' n:',(omega(i,j),j=1,nt)
            if (ispool.ne.0) then
               write (kanal0,223) ' a:',(domega(i,j,k),j=1,nt)
               write (kanal0,223) ' n:',(omega(i,j),j=1,nt)
            end if
223         format(a3,20f8.4)
224      continue
         parm(nb+k)=parm(nb+k)-eps
         if (ipause.ge.0) write (*,227) devmax
227      format(/' ==>  max absolute error = ',d15.6)
         if (ipause.ge.0) pause 'hit enter for more....'
220   continue
      return
99    write (*,*) 'error in dcomeg: inadmissible omega matrix'
      if (ispool.ne.0) write (kanal0,*)
     *            'error in dcomeg: inadmissible omega matrix'
      return
      end
c
c
      subroutine fdiff(parm,np,xlf,grad,*)
c--------------------------------------------------------------------
c     numerical first derivatives
c--------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension parm(np),grad(np)
c
ccccc (dummy statement, do not update xlf!!)
ccccc dummy=xlf
c
      call func(parm,np,xlf1,*98)
cccc             write (*,'(a,15f10.6)') 'fdiff:',(parm(j),j=1,np),xlf1
      do 100 i=1,np
         ai=parm(i)
         ei=dmax1(dabs(ai)*.1d-3,.1d-5)
12       parm(i)=ai+ei
         call func(parm,np,xlf2,*15)
cccc             write (*,'(a,15f10.6)') 'fdiff:',(parm(j),j=1,np),xlf2
         grad(i)=(xlf2-xlf1)/ei
         goto 10
15       ei=ei/10.
         if (ei.lt..1d-10) goto 99
         goto 12
10       parm(i)=ai
100   continue
      return
98    write (*,*) 'error in fdiff: first func call'
      return 1
99    write (*,*) 'error in fdiff: can''t find suffiently small step'
      return 1
      end
c
      subroutine sdiff(parm,np,grad,hess,*)
c--------------------------------------------------------------------
c     numerical second derivatives from analytic first derivatives
c--------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
      dimension parm(np),grad(np),grad1(maxnp),hess(np,np)
      character*21 tdate
      character*20 runtag
      common / dtag  / runtag
      common / doutp / kanal0,ispool,ipause
c
      idiff2=0
      call deriv(parm,np,xlf,grad1,hess,idiff2,*98)
      call tstamp(tdate)
      fpnorm=dsqrt(dotp(grad1,np,grad1))
      write (*,101) 0,fpnorm,tdate,runtag
      if (ispool.ne.0) write (kanal0,101) 0,fpnorm,tdate,runtag
101   format(' sdiff(',i2,'): fpnorm=',d12.6,'; ',a21,'; ',a20)
c
      do 100 i=1,np
         ai=parm(i)
         ei=dmax1(dabs(ai)*.1d-3,.1d-5)
12       parm(i)=ai+ei
         idiff2=0
         call deriv(parm,np,xlf,grad,hess,idiff2,*15)
         call tstamp(tdate)
         fpnorm=dsqrt(dotp(grad,np,grad))
         write (*,101) i,fpnorm,tdate,runtag
         if (ispool.ne.0) write (kanal0,101) i,fpnorm,tdate,runtag
         do 13 j=1,np
13          hess(i,j)=(grad(j)-grad1(j))/ei
         goto 10
15       ei=ei/10.
         if (ei.lt..1d-10) goto 99
         goto 12
10       parm(i)=ai
100   continue
c
      do 20 i=1,np
      do 20 j=i+1,np
         h1=hess(i,j)
         h2=hess(j,i)
         h0=0.5*h1+0.5*h2
         hess(i,j)=h0
         hess(j,i)=h0
20    continue
c
      return
98    write (*,*) 'error in sdiff: first deriv call'
      return 1
99    write (*,*) 'error in sdiff: can''t find suffiently small step'
      return 1
      end
c
c
      subroutine promeg(parm,np,kanal)
c----------------------------------------------------------------------
c     print omega 
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
c
      dimension parm(np)
      character*8 fmt
      logical lderiv
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / domeg / omega(maxnt,maxnt),domega(maxnt,maxnt,0:maxns)
      common / doutp / kanal0,ispool,ipause
c
      lderiv=.false.
      call omega1(parm,np,*99)
      fmt='(15f5.2)'
      if (nt.le.12) fmt='(12f6.3)'
      if (nt.ge.16) fmt='(20f4.1)'
      if (kanal.eq.0) then
         write (*,'(//a)') ' omega matrix='
         do 15 i=1,nt
15          write (*,fmt) (omega(i,j),j=1,nt)
         call snooze(ipause)
      else
         write (kanal,'(//a)') ' omega matrix='
         do 25 i=1,nt
25          write (kanal,fmt) (omega(i,j),j=1,nt)
      end if
      return
99    write (*,*) 'error in promeg: inadmissible omega matrix'
      if (ispool.ne.0) write (kanal0,*)
     *            'error in promeg: inadmissible omega matrix'
      return
      end
c
      subroutine prtmat(text,amat,n,m,ndim)
c--------------------------------------------------------------------
c     print matrix
c--------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      character*(*) text
      dimension amat(ndim,m)
      common / doutp / kanal0,ispool,ipause
c
      write (*,1) text,n,m,ndim
      if (ispool.gt.0) write (kanal0,1) text,n,m,ndim
1     format(' matrix ',a6,' [',i2.2,',',i2.2,'] [',i2.2,',*]')
      do 10 i=1,n
         write (*,'(9f8.3)') (amat(i,j),j=1,m)
         if (ispool.gt.0) write (kanal0,'(9f8.3)') (amat(i,j),j=1,m)
10    continue
      return
      end
c
c
      subroutine smlout(parm,hinv,np,
     *                  xlf,ier,iter,ival,ialg,fpnorm,grdir,nsimul)
c--------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
      dimension    parm(np),hinv(np,np)
      character*21 tdate
      character*20 text(3)
      character*52 text2(8)
      character*4  text3(4)
      character*8  pname(maxnpp)
      character*20 runtag
      character*40 spfile,pafile
c
      common pvar(maxnpp),lmap(maxalt)
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / dmode / method,model,modran,modar1,modmnp
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / dname / pname
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / doutp / kanal0,ispool,ipause
      common / dtag  / runtag
      common / dfile / spfile,pafile
c
      data text/'multinomial logit   ',
     *          'clark approximation ',
     *          'simulated likelihood'/
c
      data text2/
     *  'iid across periods, iid across alternatives         ',
     *  'iid across periods, correlated across alternatives  ',
     *  'random effects, iid across alternatives             ',
     *  'random effects, correlated across alternatives      ',
     *  'ar1 errors, iid across alternatives                 ',
     *  'ar1 errors, correlated across alternatives          ',
     *  'random effects + ar1 errors, iid across alternatives',
     *  'random effects + ar1 errors, correlated across alts.'/
c
      data text3/' ???',' dfp','bfgs','bhhh'/
c
      kanal8=8
      open(kanal8,file=pafile,status='old',access='append')
      if (ispool.ne.0)
     *   open(kanal0,file=spfile,status='old',access='append')
      itext=3
      if (method.eq.1) itext=1
      if (method.eq.2) itext=2
      write (*,84) model,text2(model),text(itext),text3(ialg)
      if (ispool.ne.0) write (kanal0,84) 
     *             model,text2(model),text(itext),text3(ialg)
84    format(//' ===================================================='
     *        /' sml 2.1  multinomial-multiperiod probit model    ['
     *                                                         ,i1,']'
     *        /' ===================================================='
     *        /' ',a52
     *        /' ===================================================='
     *        /' estimation by: ',a20,'  [',a4,']'
     *        /' ===================================================='
     *       //' variable  estimate   transformed  std-err  t-stat'
     *        /' --------  --------   -----------  -------  ------')
c
c *** transfer deep parameters to parmp and variances to pvar ***
c     undo transformation for sigmas
c
      npp=ns5+nalt
      j=0
      do 171 i=1,npp
         pvar(i)=0.0
         if (loc(i).ne.0) then
            j=j+1
            parmp(i)=parm(j)
            pvar(i)=-hinv(j,j)
         end if
171   continue
c
c *** output of parameters ***
c
180   do 181 i=1,npp
         m=loc(i)
         std = 0.0
         tst = 0.0
         if (m.gt.0) then
            if (pvar(i).gt.0.) then
               std = dsqrt(pvar(i))
               tst = parmp(i)/std
            else
               std = 99.99
               tst = 99.99
            end if
         end if
c
c        *** undo transformation ***
c
         if (i.le.nxy) then
c        ---beta: [divide logit estimates by sqrt(pi/6)]
            pt=parmp(i)
            if (method.eq.1) pt=0.7797*pt
         else if ((i.le.ns2).or.(i.gt.ns3.and.i.le.ns4)) then
c        ---std.dev's nu and alfa:
            pt=exp(parmp(i))
         else
c        ---rho, correl nu and alfa:
            pt=(exp(parmp(i))-1)/(exp(parmp(i))+1)
         end if
c
c        *** write ***
c
185      if (m.gt.0) write (*,186) pname(i),pt,parmp(i),std,tst
         if (m.eq.0) write (*,186) pname(i),pt,parmp(i)
186      format(1x,a8,f10.4,f14.4,f9.3,f8.2)
         if (ispool.ne.0) then
            if (m.gt.0) write (kanal0,186) pname(i),pt,parmp(i),std,tst
            if (m.eq.0) write (kanal0,186) pname(i),pt,parmp(i)
         end if
         j=loc(i)
         if (m.gt.0) write (kanal8,187) pname(i),j,pt,parmp(i),std,tst
         if (m.le.0) write (kanal8,187) pname(i),j,pt,parmp(i)
187      format(a8,i2,4f15.9)
181   continue
c
c *** likelihood at zero, mean square gradient
c
      call tstamp(tdate)
      write (*     ,192) model,modran,modar1,modmnp,
     *                   tdate,runtag,text(itext),text3(ialg),
     *                   xlf,nindiv,nper,nobs,nsimul,
     *                   iter,ival,fpnorm,grdir,ier
      if (ispool.ne.0) write (kanal0,192)
     *                   model,modran,modar1,modmnp,
     *                   tdate,runtag,text(itext),text3(ialg),
     *                   xlf,nindiv,nper,nobs,nsimul,
     *                   iter,ival,fpnorm,grdir,ier
      write (kanal8,192) model,modran,modar1,modmnp,
     *                   tdate,runtag,text(itext),text3(ialg),
     *                   xlf,nindiv,nper,nobs,nsimul,
     *                   iter,ival,fpnorm,grdir,ier
192   format(/' ==========================================='
     *       /' sml ver.2.1  multinomial-multiperiod probit'
     *       /' ==========================================='
     *       /' model=',i1,15x,'(ran=',i1,'  ar1=',i1,'  mnp=',i1,')'
     *       /' ==========================================='
     *       /' ',a21,', ',a20,
     *       /' estimation by: ',a20,'  [',a4,']'
     *       /' ==========================================='
     *       /' function value              ',f15.4
     *       /' number of individuals:      ',i15
     *       /' number of periods:          ',i15
     *       /' number of observations:     ',i15
     *       /' number of replications:     ',i15
     *       /' -------------------------------------------'
     *       /' total number of iterations: ',i15
     *       /' total number of evaluations:',i15
     *       /' norm of gradient:           ',d15.4
     *       /' gradient times direction:   ',d15.4
     *       /' return code:                ',i15
     *       /' ===========================================')
      close(kanal8)
      return
      end
c
c
      subroutine cwesml(parm,grad,hess,hinv,np,nwesml,*)
c-------------------------------------------------------------------
c
c     calculates correct covariance for wesml estimator
c
c     covm = inv(hess) * ( bhhh + aba' ) * inv(hess)
c
c     cwesml returns the corrected covariance matrix in hinv
c
c
c     nwesml=1: inv(hess) given in hinv
c     nwesml=2: recompute inv(hess) as finite differences of grad
c
c-------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
      dimension parm(np),grad(np),hess(np,np),hinv(np,np)
c
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / doutp / kanal0,ispool,ipause
      common / ddiff / idiff2
c
      common scra(maxnp,maxnp)
c
      if (nwesml.eq.1) then
c        ---save inverse hessian in "hess"
         do 10 i=1,np
         do 10 j=1,np
10          hess(i,j)=hinv(i,j)
      else if (nwesml.eq.2) then
c        ---compute inverse hessian in "hess"
         idiff2=0
         call sdiff(parm,np,grad,hess,*98)
ccc         call prtmat('sdiff=',hess,np,np,np)
         call gaussj(hess,grad,np,np,*99)
ccc         call prtmat('inv(sdiff)=',hess,np,np,np)
      else
ccc         write (*,*) 'warning: nwesml=',nwesml
         return 1
      end if
c
c     calculate bhhh matrix in "hinv"
c
      idiff2=4
      call deriv(parm,np,xlf,grad,hinv,idiff2,*98)
ccc      call prtmat('bhhh=',hinv,np,np,np)
c
c     pre- and post-multiply inv.hessian (hess) with bhhh (hinv)
c     store result in hinv
c
      ndim=maxnp
      call matp(hess,np,np  ,np,hinv,np,np,scra,ndim)
      call matp(scra,np,ndim,np,hess,np,np,hinv,np  )
ccc      call prtmat('wesml=',hinv,np,np,np)
c
      return
98    write (*,*) 'error in cwesml: invalid initial values'
      if (ispool.ne.0) write (kanal0,*)
     *            'error in cwesml: invalid initial values'
      return 1
99    write (*,*) 'error in cwesml: hessian singular'
      if (ispool.ne.0) write (kanal0,*)
     *            'error in cwesml: hessian singular'
      return 1
      end
c----------------------------------------------------------------------
c     smlopt: usual opt plus sml-specific itlog   / 1. september 1990
c----------------------------------------------------------------------
      subroutine opt(parm,grad,hess,hinv,np,
     *               fval,fpnorm,grdir,ier,
     *               mx,iter,acc,ialg)
c----------------------------------------------------------------------
c.pl100/.rm72/
c     driver and outer loop for nonlinear optimization routines
c     =========================================================
c     axel boersch-supan                 [version feb 16, 1990]
c
c     input:    parm,np .......initial parms and dimension
c               mx ............1=max, 2=min
c               iter ..........iteration limit
c               acc ...........accuracy for convergence criterion
c               ialg ...* only 2=davidon-fletcher-powell algorithm
c     **sml**           * and  4=newton-raphson with bhhh approximation
c
c     output:   parm,fval............ optimal parms and opt. func-value
c               grad ................ gradient
c               hess,hinv ........... hessian, inv.hessian (=-covmat)
c               fpnorm,grdir ........ norm of grad and grad*direc at opt
c               ier ................. return code
c                            -9=hessian singular
c                            -8=eigenvalues did not converge
c                            -7=numerical saddlepoint
c                            -6=cannot find improving step
c                            -5=too many function errors
c                            -4=*** out of memory ***
c                            -3=function error in gradient evaluation
c                            -2=initial values not admissible
c                            -1=iteration limit exceeded
c                             0=*** unknown error ***
c                             1=stepsize convergence achieved
c                             2=gradient convergence achieved
c                             3=function value convergence achieved
c                             4=gradient*direction convergence achieved
c
c----------------------------------------------------------------------
c
c     sml-specific:    c$include
c                      common / doutp /
c                      common / dopt9 /
c
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
c
      real*8 parm(np),grad(np),hess(np,np),hinv(np,np)
      character*1  prompt
      character*3  cmx(2)
      character*9  txt(6)
      character*21 tdate
      character*40 spfile,pafile
      character*45 text(14)
c
      common / dopt9 / a1(maxnp),a2(maxnp),del(maxnp),fpd1(maxnp),
     *                 spd1(maxnp,maxnp),p(maxnp,maxnp),delr(maxnp)
      common / doutp / kanal0,ispool,ipause
c
      common / dival / ival
      common / ddiff / idiff2
      common / dstop / istop
      common / dfile / spfile,pafile
c
      data cmx/'max','min'/
      data text/
     *          'hessian singular',
     *          'eigenvalues did not converge',
     *          'numerical saddlepoint',
     *          'cannot find improving step',
     *          'too many function errors',
     *          '*** out of memory ***',
     *          'function error in gradient evaluation',
     *          'initial values not admissible',
     *          'iteration limit exceeded',
     *          '*** unknown error ***',
     *          'stepsize convergence achieved',
     *          'gradient convergence achieved',
     *          'function value convergence achieved',
     *          'gradient*direction convergence achieved'/
      data txt/'steepdesc','davflepow','newt-raph','newt-bhhh',
     *          'quadh-1up','quadh-2up'/
c
c *** initializations ***
c
      iterl=iter
      ispd=ialg
      ialg0=99
      istop=1
      ival=0
      iteral=0
      nevalf=0
      do 10 i=1,np
10       hinv(i,i)=0.0
c
c *** outer iteration ***
c
      write (*,901) cmx(mx),np
      if (ispool.ne.0) write (kanal0,901) cmx(mx),np
901   format(/
     */' ============================================================='
     */' general nonlinear optimization routine  [version feb 16,1989]'
     */'         ***** ',a3,'imization over',i3,' parameters *****'
     */' =============================================================')
      call func(parm,np,fval,*999)
      write (*,902) fval
      if (ispool.ne.0) write (kanal0,902) fval
902   format(' initial function value =',f20.8)
900   write (*,903) txt(ialg),iterl,acc,ialg,istop
      if (ispool.ne.0) write (kanal0,903)
     *              txt(ialg),iterl,acc,ialg,istop
903   format(
     * ' -------------------------------------------------------------'
     */' default iteration parameters:',21x,'[',a9,']'
     */' iterations=',i3,',  accuracy=',f8.6,',  algorithm=',i1,
     *                                               ',  stopping=',i1
     */' -------------------------------------------------------------'
     */' algorithms:'
     */'   1 = steepdesc = method of steepest descent'
     */'   2 = davflepow = davidon-fletcher-powell algorithm'
     */'   4 = newt-bhhh = newton-raphson with bhhh approximation'
     */' -------------------------------------------------------------')
      if (ipause.lt.0) goto 905
      write (*,*)
     * ' hit return to accept default parameters [q=quit, c=change] >'
      read (*,'(a1)') prompt
      if (prompt.eq.'q'.or.prompt.eq.'q') goto 980
      if (prompt.eq.'a'.or.prompt.eq.'a') goto 905
      if (prompt.eq.' '.or.prompt.eq.'0') goto 905
      write (*,904)
904   format(' enter new values for'
     *      /'  (1) iteration limit,'
     *      /'  (2) accuracy,'
     *      /'  (3) algorithm [2=dfp, 4=bhhh]:'
     *      /'  (4) stopping [1=all crit, 2=grad only]:')
      read (*,*) iterl,acc,ialg,istop
905   call tstamp(tdate)
      write (*,906) txt(ialg),iterl,acc,ialg,istop,tdate
      if (ispool.ne.0) write (kanal0,906)
     *              txt(ialg),iterl,acc,ialg,istop,tdate
906   format(
     * ' -------------------------------------------------------------'
     */' actual iteration parameters:',22x,'[',a9,']'
     */' iterations=',i3,',  accuracy=',f8.6,',  algorithm=',i1,
     *                                               ',  stopping=',i1
     */' -------------------------------------------------------------'
     */' time = ',a21)
      iterc=0
      ival=0
      ier=0
      if (ispool.ne.0) then
         close(kanal0)
         open(kanal0,file=spfile,status='old',access='append')
      end if
      if (iterl.eq.0) goto 980
      if (ialg.le.2) then
         idiff2=2
c        ---do not initialize dfp-hessian if restart of same algorithm
         if (ialg.eq.ialg0) ier=99
         call dfpx(parm,grad,hess,hinv,np,
     *          fval,fpnorm,grdir,
     *          mx,iterl,acc,ialg,iteral,ier,
     *          a1,a2,del,fpd1,delr)
      else 
         idiff2=4
         call newton(parm,grad,hess,hinv,np,
     *            fval,fpnorm,grdir,
     *            mx,iterl,acc,ialg,iteral,ier,
     *            a1,a2,del,fpd1,delr)
      end if
      if (ier.gt.4.or.ier.lt.-9) ier=0
910   if (ier.gt.0) then
         write (*,911) text(ier+10)
         if (ispool.ne.0) write (kanal0,911) text(ier+10)
911      format(//' optimum reached: ',a45//)
      else
         write (*,912) text(ier+10)
         if (ispool.ne.0) write (kanal0,912) text(ier+10)
912      format(//' failure to reach optimum: ',a45//)
      end if
c
c *** end of outer iteration ***
c
950   iteral=iteral+iterc
      nevalf=nevalf+ival
      ialg0=ialg
      call itlog(parm,np,fval,fpnorm,grdir,iteral,'stop')
      write (*,951) ier
951   format(//' *** optimization terminated *** ier=',i3)
      if (ipause.lt.0) goto 980
      write (*,*) 'hit return to quit, no to continue >'
      read (*,'(a1)') prompt
      if (prompt.eq.' '.or.prompt.eq.'0') goto 980
      goto 900
999   ier=-2
      goto 910
980   ival=nevalf
      iter=iteral
      return
      end
c
      subroutine dfpx(parm,grad,hess,hinv,np,
     *                fval,fpnorm,grdir,
     *                mx,iterl,acc,ialg,iterc,ier,
     *                dir,grad1,parm1,gam,hgam)
c----------------------------------------------------------------------
c
c     nonlinear optimization routines 1 and 2   version  feb 1990
c     ===========================================================
c
c     ialg: (1) method of steepest descent
c           (2) davidon-fletcher-powell algorithm
c     mx:   maximization (1) or minimization (2)
c     acc:  accuracy for convergence criterion
c     iterl iteration limit
c
c----------------------------------------------------------------------
      implicit real*8(a-h,o-z)
      implicit integer (i-n)
      real*8 parm(np),grad(np),hess(np,np),hinv(np,np),
     *       dir(np),parm1(np),grad1(np),gam(np),hgam(np)
      character*4 text
c
      common / doutp / kanal0,ispool,ipause
      common / dstop / istop
c
      data step1,stpmin,stracc,rsmall,absmal
     *     /1.d0,1.d-10,1.0d-5,1.d-14,1.d-35/
c
c     start: initialize, initial function value
c            (unless restart with dfp)
c
      if (ier.eq.99) then
         ier=0
         text='more'
         goto 2000
      end if
c
      pm1=float(3-2*mx)
      ier=0
      text='strt'
      st0=step1
      grdir=0.0
      delf=2*fval
      f1=-pm1*1.d20
      f2=-pm1*1.d20
      f3=-pm1*1.d20
      iterl=iterc+iterl
      if (iterc.gt.0) goto 2000
c
c     initial grad and hess evaluation
c
      ialg1=1
      call deriv(parm,np,fval,grad,hess,ialg1,*9018)
c
c     initial h-matrix for dfp and steepest descent
c
      do 20 i=1,np
         do 21 j=1,np
21          hinv(i,j)=0.
         hinv(i,i)=-pm1
20    continue
c
c     main iteration loop
c
2000  if (iterc.ge.iterl) goto 9015
      gradn=dotp(grad,np,grad)
      fpnorm=dsqrt(gradn)
      call itlog(parm,np,fval,fpnorm,grdir,iterc,text)
      iterc=iterc+1
c
c     check convergence by norm of gradient
c
      if (fpnorm.le.acc) goto 9002
c
c     check convergence by function value (stationarity in three steps)
c
      fmove=dabs(fval-f1)+dabs(f1-f2)+dabs(f2-f3)
      if (fmove.lt.fval*acc.and.istop.ne.2) goto 9003
      f3=f2
      f2=f1
      f1=fval
c
c     conjugate gradient direction
c
      if (ialg.eq.1) goto 110
      do 105 i=1,np
         sum=0.d0
         do 106 j=1,np
106         sum=sum+hinv(i,j)*grad(j)
105      dir(i)=-sum
c
c     watch for singularity
c
109   grdir=dotp(dir,np,grad)
      if (grdir*pm1.gt.stracc**2) goto 120
c
c     steepest descent direction
c
110   grdir=gradn*pm1
      text='grad'
      do 115 i=1,np
115      dir(i)=pm1*grad(i)
c
c     save old parameter and function values
c
120   fval1=fval
      grdir1=grdir
      call move(parm,np,parm1)
      call move(grad,np,grad1)
c
c     linesearch
c
      z=1.0
      if (iterc.le.np) z=dabs(st0)
      st0=dmax1(stpmin,dmin1(z,dabs(2.d0*delf/grdir)))
      call linesr(parm,np,fval,dir,pm1,st0,acc,gam,*9016)
c
c     check convergence by parameter movement
c
130   ier=1
      do 131 i=1,np
         if (dabs(parm(i)-parm1(i)).gt.acc*dabs(parm(i))+rsmall) ier=0
131   continue
      if (istop.eq.2) ier=0
c
c     accept step
c
140   delf=2.d0*(fval-fval1)
145   call deriv(parm,np,fval,grad,hess,ialg,*9018)
      grdir=dotp(dir,np,grad)
      if (ialg.eq.1) goto 2000
c
c     davidon-fletcher-powell
c
170   z=-1.d0
      call step(grad,np,z,grad1,gam)
      dg=st0*(grdir-grdir1)
      if (dg*pm1.ge.-absmal) goto 190
c
c     --update h-matrix
c
180   do 181 i=1,np
         sum=0.d0
         do 182 j=1,np
182         sum=sum+hinv(i,j)*gam(j)
181      hgam(i)=sum
      ghg=dotp(gam,np,hgam)
      if (ghg*pm1.ge.-absmal) ghg=.01*dg
      ca= 1.d0/dg
      cb=-1.d0/ghg
      c4= 1.d0/(grdir1-dg/st0**2)
      c3=-dg*c4/st0
      call step(gam,np,c3,grad1,parm1)
      c3= st0/(c3*dg)
      do 183 i=1,np
183      dir(i)=st0*dir(i)
      if (dg*pm1.gt.ghg*pm1) goto 185
c
c     --broyden variant of dfp formula
c
      text='dfpb'
      ca=1.0/(dg+ghg)
      cb=-ca
      call move(gam,np,parm1)
      c3=1.d0/dg
      c4=1.d0/grdir1
      z=1.0+ghg/dg
      do 184 i=1,np
184      dir(i)=dir(i)*z-hgam(i)
c
c     --original dfp formula
c
185   text='dfph'
      do 186 i=1,np
         diri=dir(i)*ca
         hgami=cb*hgam(i)
         do 186 j=i,np
            hinv(i,j)=hinv(i,j)+diri*dir(j)+hgami*hgam(j)
            hinv(j,i)=hinv(i,j)
186      continue
      if (pm1*grdir.ge.0.5*grdir1*pm1) st0=2.0*st0
      if (ier.eq.1) goto 9001
      goto 2000
c
c     --use old h-matrix
c
190   text='dfpo'
      grdir1=grdir
      st0=4.d0*st0
      if (ier.eq.1) goto 9001
      goto 2000
c
c     successful exits
c
9001  ier=1
      goto 9000
9002  ier=2
      goto 9000
9003  ier=3
      goto 9000
c
c     failure exits
c
9015  ier=-1
      goto 9000
9016  ier=-6
      goto 9000
9018  ier=-3
      goto 9000
9000  return
      end
c
      subroutine newton(parm,grad,hess,hinv,np,
     *                  fval,fpnorm,grdir,
     *                  mx,iterl,acc,ialg,iterc,ier,
     *                  dir,grad1,parm1,gam,hgam)
c----------------------------------------------------------------------
c
c     nonlinear optimization routines 3 and 4   version  feb 1990
c     ===========================================================
c
c     ialg: (3) newton-raphson
c           (4) bhhh
c     mx:   maximization (1) or minimization (2)
c     acc:  accuracy for convergence criterion
c     iterl iteration limit
c
c----------------------------------------------------------------------
      implicit real*8(a-h,o-z)
      implicit integer (i-n)
      real*8 parm(np),grad(np),hess(np,np),hinv(np,np),
     *       dir(np),parm1(np),grad1(np),gam(np),hgam(np)
      character*4 text,meth(4)
c
      common / doutp / kanal0,ispool,ipause
      common / dstop / istop
c
      data meth/'grad','dfpx','newt','bhhh'/
      data step1,stpmin,rsmall
     *     /1.d0,1.d-10,1.d-14/
c
c     start: initialize, initial function value
c
      hgam(1)=1.d0
      pm1=float(3-2*mx)
      ier=0
      text='strt'
      st0=step1
      grdir=0.0
      delf=2*fval
      f1=-pm1*1.d20
      f2=-pm1*1.d20
      f3=-pm1*1.d20
      iterl=iterc+iterl
      if (iterc.gt.0) goto 2000
c
c     main iteration loop
c
2000  if (iterc.ge.iterl) goto 9015
c
c     grad and hess evaluation
c
      text=meth(ialg)
      call deriv(parm,np,fval,grad,hess,ialg,*9018)
c
c     newton-raphson direction: solve dir = hess(-1)*grad
c
      do 151 i=1,np
         dir(i)=grad(i)
         do 151 j=1,np
151         hinv(i,j)=hess(i,j)
      call gaussj(hinv,dir,np,np,*9019)
      grdir=dotp(dir,np,grad)
c
c     iteration log
c
      gradn=dotp(grad,np,grad)
      fpnorm=dsqrt(gradn)
      call itlog(parm,np,fval,fpnorm,grdir,iterc,text)
      iterc=iterc+1
c
c     check convergence by norm of gradient and by gradient*direction
c
      if (fpnorm.le.acc.and.istop.ne.2) goto 9002
      if (dabs(grdir).le.acc.and.istop.ne.2) goto 9004
c
c     check convergence by function value (stationarity in three steps)
c
      fmove=dabs(fval-f1)+dabs(f1-f2)+dabs(f2-f3)
      if (fmove.lt.fval*acc.and.istop.ne.2) goto 9003
      f3=f2
      f2=f1
      f1=fval
c
c     save old parameter and function values
c
120   fval1=fval
      grdir1=grdir
      call move(parm,np,parm1)
      call move(grad,np,grad1)
c
c     linesearch
c
      z=1.0
      if (iterc.le.np) z=dabs(st0)
      st0=dmax1(stpmin,dmin1(z,dabs(2.d0*delf/grdir)))
      call linesr(parm,np,fval,dir,pm1,st0,acc,gam,*9016)
c
c     check convergence by parameter movement
c
130   ier=1
      do 131 i=1,np
         if (dabs(parm(i)-parm1(i)).gt.acc*dabs(parm(i))+rsmall) ier=0
131   continue
      if (istop.eq.2) ier=0
      if (ier.eq.1) goto 9001
c
c     next iteration
c
140   delf=2.d0*(fval-fval1)
      goto 2000
c
c     successful exits
c
9001  ier=1
      goto 9000
9002  ier=2
      goto 9000
9003  ier=3
      goto 9000
9004  ier=4
      goto 9000
c
c     failure exits
c
9015  ier=-1
      goto 9000
9016  ier=-6
      goto 9000
9018  ier=-3
      goto 9000
9019  ier=-9
      goto 9000
9000  return
      end
c
      subroutine linesr(parm,np,fval,dir,pm1,st0,acc,parm1,*)
c----------------------------------------------------------------------
c     linesearch:  computes new parm and fval along dir.
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension parm(np),dir(np),parm1(np)
      call brack(parm,np,fval,dir,pm1,st0,st1,f1,st2,f2,parm1,iedge)
      if (iedge.eq.0) then
         call quadin(parm,np,fval,dir,pm1,st0,st1,f1,st2,f2,acc,parm1)
      else
         call move(parm1,np,parm)
      end if
      return
      end
c
      subroutine brack(parm,np,f0,dir,pm1,st0,st1,f1,st2,f2,parm1,iedge)
c----------------------------------------------------------------------
c     linesearch: find initial bracket [st0,st1,st2]
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension parm(np),dir(np),parm1(np)
c
c     compute second point: take initial stepsize in search direction
c
      st1=st0
      st0=0.
10    call step(parm,np,st1,dir,parm1)
      call func(parm1,np,f1,*12)
      goto 20
c     ---failure of function evaluation
12       st1 = st1/2.d0
         if (st1.lt.1.d-31) goto 99
         goto 10
c
c     if new point is better, swap
c
20    iedge = 0
      if (f1*pm1.gt.f0*pm1) then
         f=f1
         f1=f0
         f0=f
         st=st1
         st1=st0
         st0=st
      end if
c
c     compute third point: deepen step in improving direction
c
      st2=2.0*st0-st1
30    call step(parm,np,st2,dir,parm1)
      call func(parm1,np,f2,*32)
cc      write (*,*) 'brack: fval=',f2
      if (f2*pm1.le.f0*pm1) goto 90
      goto 40
c     ---failure of function evaluation
32       shold = st2
         st2=(st2+st0)/2.
         if (shold.eq.st2) goto 95
         goto 30
c
c     if improvement, deepen bracket even more
c
40    f1=f0
      f0=f2
      st1=st0
      st0=st2
      st2=3.*st0-2.*st1
      goto 30
c
c     iedge indicates that the linesearch lead to points where the
c     function was undefined (do not use inv. quadratic interpolation)
95    iedge = 1
c
c     normal exit: no improvement found, use inv. quadr. interpolation
90    call step(parm,np,st0,dir,parm1)
      return
c
c     failure
99    iedge = 1
      return
      end
c
      subroutine quadin(parm,np,f0,dir,pm1,st0,st1,f1,st2,f2,acc,parm1)
c----------------------------------------------------------------------
c     linesearch: refine given bracket [st0,st1,st2] by inverse
c                 quadratic interpolation
c----------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      dimension parm(np),dir(np),parm1(np)
      data nlnsr/10/
c
      improv=0
      ilnsr=0
      dd=dotp(dir,np,dir)
      if (dd.eq.0) goto 101
      stpac=acc/dd
c
c     iteration (max. two improvements or nlnsr steps)
c
10    if (ilnsr.gt.nlnsr) goto 100
      ilnsr=ilnsr+1
c
c     inverse quadratic interpolation
c
      stp=st0**2*(f1-f2)+st1**2*(f2-f0)+st2**2*(f0-f1)
      std=(2.*(st0*(f1-f2)+st1*(f2-f0)+st2*(f0-f1)))
c
c     check for stepsize convergence
c
      if (dabs(std).lt.1.d-31) goto 100
      stp=stp/std
      if (dabs(stp-st0).le.stpac) goto 100
c
c     define new bracket
c
      call step(parm,np,stp,dir,parm1)
      call func(parm1,np,fnew,*100)
cc      write (*,*) 'quadi: fval=',fnew
      if (fnew*pm1.gt.f0*pm1) then
         improv=1
         if ((stp-st0)*(st1-st0).gt.0.) then
            st2=st0
            f2=f0
            st0=stp
            f0=fnew
         else
            st1=st0
            f1=f0
            st0=stp
            f0=fnew
         end if
      else
         if (improv.eq.1) goto 100
         if ((stp-st0)*(st1-st0).gt.0.) then
            st1=stp
            f1=fnew
         else
            st2=stp
            f2=fnew
         end if
      end if
      goto 10
c
c     exit
c
100   call step(parm,np,st0,dir,parm1)
      call move(parm1,np,parm)
      return
101   write (*,*) 'quadin: dir=0'
      call move(parm1,np,parm)
      return
      end
c
c
c
      subroutine itlog(parm,np,fval,fpnorm,grdir,iter,text)
c---------------------------------------------------------------------
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
ccc$include:'sml.dim'
      include 'ssmlmnp.inc'
c$include:'sml.dim'
c.......................................
c.......................................
      dimension parm(np),parm1(10)
      character*4 text
      character*20 runtag
      character*21 tdate
      character*40 spfile,pafile
      common / dimen / nalt,nper,nt,nt1,nindiv,nobs,nx,ny,nc,nwt,nb,ns,
     *                 nalt1,nalt2,naltc,nxy,npp,ndd,ndw,nda,ndp,mxt,mxp
      common / doffs / ns1,ns2,ns3,ns4,ns5,maps(maxalt,maxalt),
     *                 locs(maxnpp)
      common / dparm / parmp(0:maxnpp),loc(maxnpp),loc1(maxnp)
      common / doutp / kanal0,ispool,ipause
      common / dival / ival
      common / dtag  / runtag
      common / dfile / spfile,pafile
c
      if (text.eq.'quit') return
      call tstamp(tdate)
      write (*,10) text,iter,fval,fpnorm,grdir,ival
      if (nb.gt.0) write (*,11) (parm(i),i=1,nb)
      if (ns.gt.0) then
c        *** undo transformation ***
         do 5 k=1,ns
            k1=nb+k
            i1=loc1(k1)
            if ((i1.le.ns2).or.(i1.gt.ns3.and.i1.le.ns4)) then
c           ---std.dev's nu and alfa:
               parm1(k)=exp(parm(k1))
            else
c           ---rho, correl nu and alfa:
               parm1(k)=(exp(parm(k1))-1)/(exp(parm(k1))+1)
            end if
5        continue
         write (*,12) (parm1(k),k=1,ns)
      end if
      write (*,13) tdate,runtag
      if (ispool.ne.0) then
         write (kanal0,10) text,iter,fval,fpnorm,grdir,ival
         if (nb.gt.0) write (kanal0,11) (parm(i),i=1,nb)
         if (ns.gt.0) write (kanal0,12) (parm1(k),k=1,ns)
         write (kanal0,13) tdate,runtag
         close(kanal0)
         open(kanal0,file=spfile,status='old',access='append')
      end if
10    format(/1x,a4,'-iteration',i3,':',
     *        ' lik=',d12.5,', gradnm=',d9.3,', grad*dir=',d9.3,i5)
11    format(' beta =',d11.4,5d12.4:/10(6x,6d12.4:/))
12    format(' sigma=',d11.4,5d12.4:/10(6x,6d12.4:/))
13    format(' time = ',a21,', ',a20)
      return
      end

c----------------------------------------------------------------------c
c     system-dependent routines                                        c
c     If you cannot find the correct calls for your system,            c
c     use the generic ones, currently commented out by CDUMMY          c
c                                                                      c
c----------------------------------------------------------------------c
      subroutine getdat(iy,im,id)
c     system-dependent routine: SUN/UNIX version plus VAH library
      integer*2 iy,im,id
      integer*4 itim(5),idat(3)
      call timdat(itim,idat)
      id=idat(1)
      im=idat(2)
      iy=idat(3)
cvah  since timdat/unix returns year in two digits
      iy=iy+1900
      return
      end

      subroutine gettim(ih,im,is,ihs)
c     system-dependent routine: SUN/UNIX version plus VAH library
      integer*2 ih,im,is,ihs
      integer*4 itim(5),idat(3)
      call timdat(itim,idat)
      ih=itim(1)
      im=itim(2)
      is=itim(3)
      ihs=itim(4)
      return
      end

cdummy       subroutine gettim(ih,im,is,ihs)
cdummy c     returns 0 for hour, minute, and 10th of sec, and adds 1 to
cdummy c     previous value of second.
cdummy c     ideally, should return the time as
cdummy c     ih=hour, im=minute, is=second, ihs=hundredth of a sec
cdummy       integer*2 ih,im,is,ihs
cdummy       integer itemp
cdummy       save itemp
cdummy       ih=0
cdummy       im=0
cdummy       itemp=itemp+1
cdummy       is=itemp
cdummy       ihs=0
cdummy       return
cdummy       end
cdummy       subroutine getdat(iy,im,id)
cdummy       integer*2 iy,im,id
cdummy c     ideally, should return date as
cdummy c     iy=year, im=month, id=day of the month
cdummy       iy=1991
cdummy       im=1
cdummy       id=1
cdummy       return
cdummy       end
cdummy       subroutine parsecl(cline,istrt,iend,narg)
cdummy c     ideally, returns the command line when the program was
cdummy c     called as cline.  It also parses cline and returns
cdummy c     number of distinct arguments, narg, as well as integer
cdummy c     arrays of length narg, specifying where each argument
cdummy c     begins (istrt) and ends (iend)
cdummy       character*(*) cline
cdummy       integer*4 istrt,iend,narg
cdummy       istrt=1
cdummy       iend=1
cdummy       narg=0
cdummy       return
cdummy       end
