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