C devivatives of pij=prob{Y(t)=1|y(t-2)=i,y(t-1)=j}
C wr to theta(t),theta(t-1), theta(t-2), psi1, psi2
C j=y(i0) (j=y_(t-2)), j1=y(i1)  (j1=y_(t-1))

      subroutine deriv2(theta,psi1,psi2,n,t,j,j1,der)
      implicit double precision (a-h,o-z)
      integer t,j,j1,j2,j3,j4
      dimension theta(n),der(5)
      th = theta(t)
      th1 = theta(t-1) 
      th2 = theta(t-2)
      ps1 = psi1-1 
      ps2 = psi2-1 

  
      if (dabs(ps1).gt.1.0e-6.and.dabs(ps2).gt.1.0e-6) then

      d0=dsqrt(1+ps1*(psi1*(th-th1)**2-(th+th1)**2+2*(th+th1)))
      d0th=(ps1*(psi1*(th-th1)-(th+th1)+1))/d0
      d0th1=(ps1*(-psi1*(th-th1)-(th+th1)+1))/d0
      d0th2=0
      d0psi1=((2*psi1-1)*(th-th1)**2-(th+th1)**2+2*(th+th1))/(2*d0)
      d0psi2=0   
      d=dsqrt(1+ps1*(psi1*(th1-th2)**2-(th1+th2)**2+ 2*(th1+th2)))
      dth=0
      dth1=ps1*(psi1*(th1-th2)-(th1+th2)+1)/d
      dth2=ps1*(-psi1*(th1-th2)-(th1+th2)+1)/d
      dpsi1=((2*psi1-1)*(th1-th2)**2-(th1+th2)**2+2*(th1+th2))/(2*d)
      dpsi2=0

      d1=dsqrt(th1**2+(ps2/(4*ps1**2))*(ps2*(d0-d+ps1*(th2-th))**2-
     *4*(d-1)*(d0-1)+ 4*ps1*((d-1)*th+(d0-1)*th2)+
     * 4*ps1**2*(th1**2-th*th2)))
      d1th=(1/(2*d1))*(1/(4*ps1**2))*(2*ps2**2*(d0th-ps1)*
     *(d0-d+ps1*(th2-th))+4*ps1**2*ps2*(-th2)+4*ps1*ps2*
     *   (d-1+d0th*th2)-4*ps2*(d0th*(d-1)))
      d1th1=(th1/d1)+(1/(2*d1))*(1/(4*ps1**2))*(2*ps2**2*(d0th1-dth1)*
     *  (d0-d+ps1*(th2-th))+ 4*ps1**2*ps2*(2*th1)+4*ps1*ps2*
     *  (dth1*th+d0th1*th2)-4*ps2*(dth1*(d0-1)+d0th1*(d-1)))
      d1th2=(1/(2*d1))*(1/(4*ps1**2))*(2*ps2**2*(-dth2+ps1)*
     *   (d0-d+ps1*(th2-th)) +4*ps1**2*ps2*(-th)+4*ps1*ps2*
     * (dth2*th+d0-1)-4*ps2*(dth2*(d0-1)))
      d1psi1=(1/(2*d1))*(1/(4*ps1**2))*(2*ps2**2*(d0psi1-dpsi1+th2-th)*
     *  (d0-d+ps1*(th2-th))+8*ps1*ps2*(th1**2-th*th2)+4*ps2*
     * ((d-1)*th+(d0-1)*th2)+4*ps1*ps2*(dpsi1*th+d0psi1*th2)-4*ps2*
     * (dpsi1*(d0-1)+d0psi1*(d-1)))-(1/(2*d1))*(d1**2-th1**2)*(2/ps1)
      d1psi2=(1/(2*d1))*(1/(4*ps1**2))*(2*ps2*(d0-d+ps1*(th2-th))**2+
     * 4*ps1**2*(th1**2-th*th2)+4*ps1*((d-1)*th+(d0-1)*th2)-
     * 4*(d-1)*(d0-1))
      
      d2=dsqrt((1-th1)**2+(ps2/(4*ps1**2))*(ps2*(d-d0+ps1*(th2-th))**2- 
     *4*(1-d)*(1-d0)+ 4*ps1*((d-1)*(1-th)+(d0-1)*(1-th2))+
     *  4*ps1**2*((1-th1)**2+(1-th)*(th2-1)) ) )
      d2th=(1/(2*d2))*(1/(4*ps1**2))*(2*ps2**2*(-d0th-ps1)*
     *(d-d0+ps1*(th2-th))+ 4*ps1**2*ps2*(1-th2)+
     *4*ps1*ps2*(1-d+d0th*(1-th2))- 4*ps2*(d0th*(d-1)) )      
      d2th1=((th1-1)/d2)+(1/(2*d2))*(1/(4*ps1**2))*(2*ps2**2*
     *  (dth1-d0th1)*(d-d0+ps1*(th2-th))+4*ps1**2*ps2*(2*(th1-1))+ 
     *  4*ps1*ps2*(dth1*(1-th)+d0th1*(1-th2))-
     * 4*ps2*(dth1*(d0-1)+d0th1*(d-1)))
      d2th2=(1/(2*d2))*(1/(4*ps1**2))*(2*ps2**2*(dth2+ps1)*
     *(d-d0+ps1*(th2-th))+ 4*ps1**2*ps2*(1-th)+4*ps1*ps2*
     *(dth2*(1-th)+1-d0)- 4*ps2*(dth2*(d0-1)) )
      d2psi1=(1/(2*d2))*(1/(4*ps1**2))*(2*ps2**2*(dpsi1-d0psi1+th2-th)*
     *   (d-d0+ps1*(th2-th))+8*ps1*ps2*((th1-1)**2+ (1-th)*(th2-1))+
     *   4*ps2*((d-1)*(1-th)+(d0-1)*(1-th2))+4*ps1*ps2*
     *(dpsi1*(1-th)+d0psi1*(1-th2))-   4*ps2*(dpsi1*(d0-1)+
     * d0psi1*(d-1)) )-(1/(2*d2))*(d2**2-(1-th1)**2)*(2/ps1)
      d2psi2=(1/(2*d2))*(1/(4*ps1**2))*(2*ps2*(d-d0+ps1*(th2-th))**2+
     * 4*ps1**2*((th1-1)**2+(1-th)*(th2-1))+ 4*ps1*((d-1)*(1-th)+
     * (d0-1)*(1-th2))-4*(1-d)*(1-d0))
C80      
      j2 = 2*j1-1
      j3 = 2*j-1
      j4 = 2*(j+j1-2*j*j1)-1
      A=2*ps2*(j4*(d-1)+ps1*(2-2*(j+j1-j*j1)+j2*th1+j3*th2))
      B=ps2*(2*j*j2-j2*d0+j4*d+ps1*(th+2*j*j2*th1+j3*th2))-
     *    2*j1*j3*ps1*(d1-th1)+2*(j1-1)*j3*ps1*(th1-1+d2)
       dpth=(1/A)*(ps2*(-j2*d0th+ps1)-2*j1*j3*ps1*d1th+
     * 2*(j1-1)*j3*ps1*d2th)
       dpth1=(1/A)*(ps2*(-j2*d0th1+j4*dth1+2*j*j2*ps1)-
     * 2*j1*j3*ps1*(d1th1-1)+2*(j1-1)*j3*ps1*(1+d2th1))- 
     *(B/A**2)*(2*ps2*(j4*dth1+j2*ps1))
       dpth2=(1/A)*(ps2*(j4*dth2+j3*ps1)-2*j1*j3*ps1*d1th2+
     * 2*(j1-1)*j3*ps1*d2th2)-(B/A**2)*(2*ps2*(j4*dth2+j3*ps1))      
      dppsi1=(1/A)*(ps2*(-j2*d0psi1+j4*dpsi1+(th+2*j*j2*th1+j3*th2))-
     *  2*j1*j3*(d1-th1)-2*j1*j3*ps1*d1psi1+2*(j1-1)*j3*(th1-1+d2)+
     *  2*(j1-1)*j3*ps1*d2psi1)-
     *  (B/A**2)*2*ps2*(j4*dpsi1+(2-2*(j+j1-j*j1)+j2*th1+j3*th2))
      dppsi2=(1/A)*(2*j*j2-j2*d0+j4*d+ps1*(th+2*j*j2*th1+j3*th2)-
     *  2*j1*j3*ps1*d1psi2+2*(j1-1)*j3*ps1*d2psi2)-
     *  (B/A**2)*2*(j4*(d-1)+ps1*(2-2*(j+j1-j*j1)+j2*th1+j3*th2))


      else if (dabs(ps1).gt.1.0e-6.and.dabs(ps2).lt.1.0e-6) then
     
      d0=dsqrt(1+ps1*(psi1*(th-th1)**2-(th+th1)**2+2*(th+th1)))
      d0th=(ps1*(psi1*(th-th1)-(th+th1)+1))/d0
      d0th1=(ps1*(-psi1*(th-th1)-(th+th1)+1))/d0
      d0th2=0
      d0psi1=((2*psi1-1)*(th-th1)**2-(th+th1)**2+2*(th+th1))/(2*d0)
      d0psi2=0   
      d=dsqrt(1+ps1*(psi1*(th1-th2)**2-(th1+th2)**2+ 2*(th1+th2)))
      dth=0
      dth1=ps1*(psi1*(th1-th2)-(th1+th2)+1)/d
      dth2=ps1*(-psi1*(th1-th2)-(th1+th2)+1)/d
      dpsi1=((2*psi1-1)*(th1-th2)**2-(th1+th2)**2+2*(th1+th2))/(2*d)
      dpsi2=0

      j2 = 2*j1-1
      j3 = 2*j-1
      j4 = 2*(j+j1-2*j*j1)-1
      A= 2*ps1*(1-j1+j2*th1)
      B=j2*(1-d0+ps1*th1)+ps1*th
C  d2d1=second derivative of delta1 in respect to psi2
      d2d1=(1/(512*ps1**6*th1**3))*(128*ps1**4*th1**2*
     *(d0-d+ps1*(th2-th))**2 -8*ps1**2*(4*ps1**2*(th1**2-
     *th*th2)+4*ps1*((d-1)*th+(d0-1)*th2)-
     *4*(1-d)*(1-d0))**2)
C  d2d2=second derivative of delta2 in respect to psi2
      d2d2=(1/(512*ps1**6*(1-th1)**3))*(128*ps1**4*(1-th1)**2*
     *(d-d0+ps1*(th2-th))**2 -8*ps1**2*(4*ps1**2*((th1-1)**2+
     *(1-th)*(th2-1))+4*ps1*((d-1)*(1-th)+(d0-1)*(1-th2))-
     *4*(1-d)*(1-d0))**2)
      dpth=(1/A)*(-j2*d0th+ps1)
      dpth1=(1/A)*(j2)*(-d0th1+ps1)-(B/A**2)*(j2*2*ps1)
      dpth2=0
      dppsi1=(1/A)*(j2*(-d0psi1+th1)+th)-
     *(B/A**2)*(2*(1-j1+j2*th1))
      dppsi2=(-j1*j3*2*ps1*d2d1+(j1-1)*j3*2*ps1*d2d2)/
     * (4*(j4*(d-1)+ps1*(2-2*(j+j1-j*j1)+j2*th1+j3*th2)))
C140

      else if (dabs(ps1).lt.1.0e-6.and.dabs(ps2).gt.1.0e-6) then
   
      d3=dsqrt(1+ps2**2*(th2-th)**2+ps2*(2*th+2*th2-4*th*th2))
      d3th=(-ps2**2*(th2-th)+ps2*(1-2*th2))/d3
      d3th1=0
      d3th2=(ps2**2*(th2-th)+ps2*(1-2*th))/d3
      d3psi1=0
      d3psi2=(ps2*(th2-th)**2+(th+th2-2*th*th2))/d3 

     
C  d1d=first derivative of delta in respect to psi1;  d1d0=first derivative of delta0 in respect to psi1
C  d2d=second derivative of delta in respect to psi1; d2d0=second derivative of delta0 in respect to psi1
C  d1d1=first derivative of delta1 in respect to psi1;  d1d2=first derivative of delta2 in respect to psi1

      d1d= -2*th1*th2+th1+th2
      d1d0=-2*th*th1+th+th1
      d2d=4*th1*th2*(1-th2)*(th1-1)
      d2d0=4*th*th1*(1-th1)*(th-1)

      d1=dsqrt(th1**2+(2*ps2**2*(d1d0-d1d+th2-th)**2+
     *    8*ps2*(th1**2-th*th2)+8*ps2*(d1d*th+d1d0*th2)-
     *    8*ps2*d1d*d1d0)/8)

      d2=dsqrt((1-th1)**2+(2*ps2**2*(d1d-d1d0+th2-th)**2+
     *    8*ps2*((1-th1)**2+(1-th)*(th2-1))+
     *    8*ps2*(d1d*(1-th)+d1d0*(1-th2))-8*ps2*d1d*d1d0)/8) 

      d1d1= (6*ps2**2*(d2d0-d2d)*(d1d0-d1d+th2-th)+
     *    12*ps2*(d2d*th+d2d0*th2)-
     *    12*ps2*(d2d*d1d0+d2d0*d1d) )/(48*d1)

      d1d2=(6*ps2**2*(d2d-d2d0)*(d1d-d1d0+th2-th)+
     *    12*ps2*(d2d*(1-th)+d2d0*(1-th2))-
     *    12*ps2*(d2d*d1d0+d2d0*d1d))/(48*d2)

      j2 = 2*j1-1
      j3 = 2*j-1
      j4 = 2*(j+j1-2*j*j1)-1
      A=2*ps2*(1-j+j3*th2)
      B=j3*(1-d3+ps2*th2)+ps2*th

      dpth=(1/A)*(-j3*d3th+ps2)
      dpth1=0
      dpth2=(1/A)*(j3)*(-d3th2+ps2)-(B/A**2)*(j3*2*ps2)
      dppsi1=0 
 
      dppsi1=(ps2*(-j2*d2d0+j4*d2d)-4*j1*j3*d1d1+
     * 4*(j1-1)*j3*d1d2)/
     * (4*ps2*(j4*d1d+2-2*(j+j1-j*j1)+j2*th1+j3*th2))-
     * ((ps2*(-j2*d1d0+j4*d1d+th+2*j*j2*th1+j3*th2)-
     * 2*j1*j3*(d1-th1)+2*(j1-1)*j3*(d2+th1-1))*2*j4*ps2*d2d)/
     * (2*(2*ps2*(j4*d1d+2-2*(j+j1-j*j1)+j2*th1+j3*th2))**2)
	
      dppsi2=(1/A)*(j3*(-d3psi2+th2)+th)-
     *(B/A**2)*(2*(1-j+j3*th2))

      else
C      if (ps1.eq.0.0.and.ps2.eq.0.0) then
C  # when psi(1)=1 and psi(2)=1      

      dpth=1
      dpth1=0
      dpth2=0
      dppsi1=(th1-j1)*(th**2-th)
      dppsi2=(th2-j)*(th**2-th)

      end if

       der(1)=dpth
       der(2)=dpth1
       der(3)=dpth2
       der(4)=dppsi1
       der(5)=dppsi2
      return
      end