C # Calculates the log-likelihood for the second
C# order dependence model

      subroutine mlik2i(logL,pij,npar,n)
      implicit double precision (a-h,o-z)
      DIMENSION x1(4500,10),theta1(4500),
     *work1(4500),y1(4500),lpsi1(2),
     *beta1(10),bt1(10),pij(n),psi(2),tpr(2),tpr1(4),
     *P0(2,2), P1(2,2), P2(2,2),P3(4,4), P4(4,4), P5(4,4), 
     *P6(4,4),Pc(4,1),Pr(2,1),Pr0(4,1),Pr1(2,1),Paux(2,2),
     *Pres(2,2)

      double precision logL,lpsi,lpsi1,pij
      integer y1,i0,i1,i2,npar,n,n0,k,iaux, iaux1, iaux2,
     *m,mpar

      COMMON/param/x1,theta1,work1,
     *y1,lpsi1,beta1,bt1,m,mpar,omega1

      psi(1) = dexp(lpsi1(1))
      psi(2) = dexp(lpsi1(2))
      psi1 = psi(1)
      psi2 = psi(2)
      ps1 = psi1-1
      ps2 = psi2-1

      call mati(x1,beta1,work1,4500,10,1,n,npar)

      do 20 i=1,n
        theta1(i) = 1/(1+dexp(-work1(i)))
   20 continue

      i0 = 1  
   25   if (y1(i0).eq.(-1)) then
      i0=i0+1
      go to 25
      end if

      n0 = n
   30   if (y1(n0).eq.(-1)) then
      n0=n0-1
      go to 30
      end if

      logL = 0
      p  = theta1(i0)
      pij(i0)=p
      logL = y1(i0)*dlog(p/(1-p))+dlog(1-p)
      if (i0.eq.n0) return
      i = i0+1
      i1 = i
   40   if (y1(i1).eq.(-1)) then
      i1=i1+1
      go to 40
      end if

C  i0 is the most recent (past) observation time
C  i1 is the next observation time
      if (i1.eq.i) then
      th1 = theta1(i)
      th2 = theta1(i-1)
      call mcpj(th1,th2,psi1,tpr)
      else 
C    (one or more intermediate missing datum)
      call mat2 (0.0D0,1.0D0,P0)
      do 50 k=(i0+1),i1
        th1 = theta1(k)
        th2 = theta1(k-1)
        call mcpj (th1,th2,psi1,tpr)
        call mat2 (tpr(1),tpr(2),P1)   
        call matp(P0,P1,P2,2,2,2)
        call matc(P2,P0,2,2)
   50 continue
      tpr(1)= P0(1,2)
      tpr(2)= P0(2,2)
      end if
      p=tpr(y1(i0)+1)
      pij(i1)=p
      logL = logL+y1(i1)*dlog(p/(1-p))+dlog(1-p)

      if (i1.eq.n0) return
      i=i1+1
   60 if (i.le.n0) then
      i2 = i
   70   if (y1(i2).eq.(-1)) then
      i2=i2+1
      go to 70
      end if

      if (i2.eq.i.and.i1.eq.(i0+1)) then
        th = theta1(i2)
        th1 = theta1(i1)
        th2 = theta1(i0)
      call mcpij(th,th1,th2,psi1,psi2,tpr1)
        p=tpr1(y1(i0)+2*y1(i1)+1)
        pij(i2)=p
        logL = logL+y1(i2)*dlog(p/(1-p))+dlog(1-p)
        i0=i1
        i1=i2
        i=i1+1

      else if (i1.ne.(i0+1).and.i2.eq.(i1+1)) then
C    (one intermediate missing datum between i0 and i1)
      call mat2 (0.0D0,1.0D0,P0)
      do 75 k=(i0+1),(i1-1)
        th1 = theta1(k)
        th2 = theta1(k-1)
        call mcpj (th1,th2,psi1,tpr)
        call mat2 (tpr(1),tpr(2),P1)   
        call matp(P0,P1,P2,2,2,2)
        call matc(P2,P0,2,2)
   75 continue
        th = theta1(i2)
        th1 = theta1(i1)
        th2 = theta1(i1-1)
        call mcpij (th,th1,th2,psi1,psi2,tpr1)
        Pr(1,1)=tpr1(2*y1(i1)+1)
        Pr(2,1)=tpr1(2*y1(i1)+2)
        call matp(P0,Pr,Pr1,2,2,1)
        p=Pr1(y1(i0)+1,1)
	pij(i2)=p
        logL = logL+y1(i2)*dlog(p/(1-p))+dlog(1-p)
        i0=i1
        i1=i2
        i=i1+1

      else if (i2.ne.(i1+1).and.i2.ne.n0) then
C    (one or more intermediate missing datum between i1 and i2)
      call mat4 (0.0D0,0.0D0,1.0D0,1.0D0,P3)
      call matp (P3,P3,P4,4,4,4)
      do 80 k=(i1+1),i2
        th = theta1(k)
        th1 = theta1(k-1)
        th2 = theta1(k-2)
        call mcpij (th,th1,th2,psi1,psi2,tpr1)
        call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5)
        call matp(P4,P5,P6,4,4,4)
        call matc(P6,P4,4,4)
   80 continue

        i3=i2+1
        th = theta1(i3)
        th1 = theta1(i3-1)
        th2 = theta1(i3-2)
        call mcpij (th,th1,th2,psi1,psi2,tpr1)
        call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5)
        call matp(P4,P5,P6,4,4,4)
        pstar=P6(2*y1(i0)+y1(i1)+1,2*y1(i2)+y1(i3)+1)
	pij(i2)=pstar
        logL = logL+dlog(pstar)
        i0=i2
        i1=i3
        i=i1+1

      else if (i2.ne.(i1+1).and.i2.eq.n0) then
C    (one intermediate missing datum between i1 and i2(last value))
      call mat4 (0.0D0,0.0D0,1.0D0,1.0D0,P3)
      call matp (P3,P3,P4,4,4,4)
      do 90 k=(i1+1),i2
        th = theta1(k)
        th1 = theta1(k-1)
        th2 = theta1(k-2)
        call mcpij (th,th1,th2,psi1,psi2,tpr1)
        call mat4 (tpr1(1),tpr1(3),tpr1(2),tpr1(4),P5)
        call matp(P4,P5,P6,4,4,4)
        call matc(P6,P4,4,4)
   90 continue
        Pc(1,1)=0
        Pc(2,1)=1
        Pc(3,1)=0
        Pc(4,1)=1
        call matp(P4,Pc,Pr0,4,4,1)
        p=Pr0(2*y1(i0)+y1(i1)+1,1)
	pij(i2)=p
        logL = logL+y1(i2)*dlog(p/(1-p))+dlog(1-p)
        i0=i1
        i1=i2
        i=i1+1
       end if
      go to 60
      end if
      return
      end