c ********************************************************************* c Compute the trace function c ********************************************************************* double precision function trace(z,m) implicit real*8(a-h,o-z) real*8 z(m,m) qq = 0. do 100 j=1,m 100 qq = qq + z(j,j) trace = qq return end subroutine findsig(th1,th2,th3,th4,th6,th7,m1,m3,sigsig,si, 6 GG1,GG2,GG3,GG4) implicit real*8(a-h,o-z) real*8 GG1(3,3),GG2(3,3),GG3(3,3),GG4(3,3), 6 sigsig(3,3),si(3,3) call zeros(GG1,3,3) call zeros(GG2,3,3) call zeros(GG3,3,3) call zeros(GG4,3,3) GG2(1,2) = 1 GG2(1,3) = 1 GG2(2,1) = 1 GG2(3,1) = 1 GG3(2,2) = 1 GG3(3,3) = 1 GG4(2,3) = 1 GG4(3,2) = 1 GG1(1,1) = 1 sigsig(1,1) = th3 sigsig(1,2) = th4 sigsig(1,3) = th4 sigsig(2,2) = th6 sigsig(3,3) = sigsig(2,2) sigsig(2,3) = th7 sigsig(2,1) = sigsig(1,2) sigsig(3,1) = sigsig(1,3) sigsig(3,2) = sigsig(2,3) call xinverse(sigsig,si,3,bigdet) return end c ******************************************************************* c Compute the likelihood score c ******************************************************************* subroutine fullmlde(th1,th2,th3,th4,th6,th7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,derder,n,msamp) implicit real*8(a-h,o-z) real*8 GG1(3,3),GG2(3,3),GG3(3,3),GG4(3,3), 6 sigsig(3,3),si(3,3) real*8 th1th2(3,1),derder(6,1),zdtzd(3,3),zdsum(3,1), 6 zd(3,1),zdzd(3,3),tempa(1,3),tempb(1,3),tempc(1,1), 6 tempd(3,3),tempf(3,3),tempg(3,3) call findsig(th1,th2,th3,th4,th6,th7,m1,m3,sigsig,si, 6 GG1,GG2,GG3,GG4) nnn = msamp th1th2(1,1) = th1 th1th2(2,1) = th2 th1th2(3,1) = th2 do 100 j=1,3 100 zd(j,1) = zdsum(j,1) - (dfloat(msamp) * th1th2(j,1)) call zeros(zdzd,3,3) do 200 j=1,3 do 180 k=1,3 zdzd(j,k) = zdtzd(j,k) - (th1th2(j,1) * zdsum(k,1)) 6 - (zdsum(j,1) * th1th2(k,1)) 6 + (dfloat(msamp) * th1th2(j,1) * th1th2(k,1)) 180 continue 200 continue call zeros(derder,6,1) call zeros(tempa,1,3) tempa(1,1) = 1 call zmmultm(tempa,si,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) derder(1,1) = tempc(1,1) 6 + ( (yotsum - ( (n - dfloat(msamp)) * th1)) / th3) call zeros(tempa,1,3) tempa(1,2) = 1 tempa(1,3) = 1 call zmmultm(tempa,si,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) derder(2,1) = tempc(1,1) call zmmultm(si,GG1,tempd,3,3,3) sszz = trace(tempd,3) derder(3,1) = -(dfloat(msamp) / 2.0d0) * sszz call zmmultm(tempd,si,tempf,3,3,3) call zmmultm(tempf,zdzd,tempg,3,3,3) sszz = trace(tempg,3) derder(3,1) = derder(3,1) + (.5 * sszz) aaaa = (yotysum - (2 * th1 * yotsum) + 6 ( (n - dfloat(msamp)) * th1 * th1)) aaaa = (aaaa / (2 * th3 * th3)) - 6 ( (n - dfloat(msamp)) / (2 * th3)) derder(3,1) = derder(3,1) + aaaa call zmmultm(si,GG2,tempd,3,3,3) sszz = trace(tempd,3) derder(4,1) = -(dfloat(msamp) / 2.0d0) * sszz call zmmultm(tempd,si,tempf,3,3,3) call zmmultm(tempf,zdzd,tempg,3,3,3) sszz = trace(tempg,3) derder(4,1) = derder(4,1) + (.5 * sszz) call zmmultm(si,GG3,tempd,3,3,3) sszz = trace(tempd,3) derder(5,1) = -(dfloat(msamp) / 2.0d0) * sszz call zmmultm(tempd,si,tempf,3,3,3) call zmmultm(tempf,zdzd,tempg,3,3,3) sszz = trace(tempg,3) derder(5,1) = derder(5,1) + (.5 * sszz) call zmmultm(si,GG4,tempd,3,3,3) sszz = trace(tempd,3) derder(6,1) = -(dfloat(msamp) / 2.0d0) * sszz call zmmultm(tempd,si,tempf,3,3,3) call zmmultm(tempf,zdzd,tempg,3,3,3) sszz = trace(tempg,3) derder(6,1) = derder(6,1) + (.5 * sszz) return end c ******************************************************************* c Compute the Hessian c ******************************************************************* subroutine fullmlhe(th1,th2,th3,th4,th6,th7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,hesshess,n,msamp) implicit real*8(a-h,o-z) real*8 GG1(3,3),GG2(3,3),GG3(3,3),GG4(3,3), 6 sigsig(3,3),si(3,3),CC1(3,3),CC2(3,3) real*8 th1th2(3,1),zdtzd(3,3),zdsum(3,1), 6 zd(3,1),zdzd(3,3),tempb(1,3),tempc(1,1), 6 tempf(3,3),tempg(3,3), 6 tempa1(1,3),tempa2(1,3),tempa1t(3,1),tempa2t(3,1) real*8 hesshess(6,6) call findsig(th1,th2,th3,th4,th6,th7,m1,m3,sigsig,si, 6 GG1,GG2,GG3,GG4) call zeros(hesshess,6,6) nnn = msamp th1th2(1,1) = th1 th1th2(2,1) = th2 th1th2(3,1) = th2 do 100 j=1,3 100 zd(j,1) = zdsum(j,1) - (dfloat(msamp) * th1th2(j,1)) call zeros(zdzd,3,3) do 200 j=1,3 do 180 k=1,3 zdzd(j,k) = zdtzd(j,k) - (th1th2(j,1) * zdsum(k,1)) 6 - (zdsum(j,1) * th1th2(k,1)) 6 + (dfloat(msamp) * th1th2(j,1) * th1th2(k,1)) 180 continue 200 continue call zeros(tempa1,1,3) tempa1(1,2) = 1 tempa1(1,3) = 1 call zeros(tempa2,1,3) tempa2(1,1) = 1 do 300 j=1,3 tempa1t(j,1) = tempa1(1,j) 300 tempa2t(j,1) = tempa2(1,j) call zmmultm(tempa1,si,tempb,1,3,3) call zmmultm(tempb,tempa1t,tempc,1,3,1) hesshess(1,1) = -nnn * tempc(1,1) call zmmultm(tempa2,si,tempb,1,3,3) call zmmultm(tempb,tempa1t,tempc,1,3,1) hesshess(2,1) = -nnn * tempc(1,1) call zmmultm(tempa1,si,tempb,1,3,3) call zmmultm(tempb,tempa2t,tempc,1,3,1) hesshess(1,2) = -nnn * tempc(1,1) call zmmultm(tempa2,si,tempb,1,3,3) call zmmultm(tempb,tempa2t,tempc,1,3,1) hesshess(2,2) = -nnn * tempc(1,1) call zmmultm(si,GG1,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa1,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(1,3) = -tempc(1,1) call zmmultm(si,GG2,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa1,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(1,4) = -tempc(1,1) call zmmultm(si,GG3,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa1,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(1,5) = -tempc(1,1) call zmmultm(si,GG4,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa1,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(1,6) = -tempc(1,1) call zmmultm(si,GG1,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa2,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(2,3) = -tempc(1,1) call zmmultm(si,GG2,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa2,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(2,4) = -tempc(1,1) call zmmultm(si,GG3,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa2,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(2,5) = -tempc(1,1) call zmmultm(si,GG4,tempf,3,3,3) call zmmultm(tempf,si,tempg,3,3,3) call zmmultm(tempa2,tempg,tempb,1,3,3) call zmmultm(tempb,zd,tempc,1,3,1) hesshess(2,6) = -tempc(1,1) hesshess(3,1) = hesshess(1,3) hesshess(4,1) = hesshess(1,4) hesshess(5,1) = hesshess(1,5) hesshess(6,1) = hesshess(1,6) hesshess(3,2) = hesshess(2,3) hesshess(4,2) = hesshess(2,4) hesshess(5,2) = hesshess(2,5) hesshess(6,2) = hesshess(2,6) do 1010 j=1,3 do 1005 k=1,3 CC1(j,k) = GG1(j,k) 1005 CC2(j,k) = GG1(j,k) 1010 continue hesshess(3,3) = zakezake(zdzd,CC1,CC2,si,nnn) do 1020 j=1,3 do 1015 k=1,3 CC1(j,k) = GG1(j,k) 1015 CC2(j,k) = GG2(j,k) 1020 continue hesshess(3,4) = zakezake(zdzd,CC1,CC2,si,nnn) do 1030 j=1,3 do 1025 k=1,3 CC1(j,k) = GG1(j,k) 1025 CC2(j,k) = GG3(j,k) 1030 continue hesshess(3,5) = zakezake(zdzd,CC1,CC2,si,nnn) do 1040 j=1,3 do 1035 k=1,3 CC1(j,k) = GG1(j,k) 1035 CC2(j,k) = GG4(j,k) 1040 continue hesshess(3,6) = zakezake(zdzd,CC1,CC2,si,nnn) do 1050 j=1,3 do 1045 k=1,3 CC1(j,k) = GG2(j,k) 1045 CC2(j,k) = GG2(j,k) 1050 continue hesshess(4,4) = zakezake(zdzd,CC1,CC2,si,nnn) do 1060 j=1,3 do 1055 k=1,3 CC1(j,k) = GG2(j,k) 1055 CC2(j,k) = GG3(j,k) 1060 continue hesshess(4,5) = zakezake(zdzd,CC1,CC2,si,nnn) do 1070 j=1,3 do 1065 k=1,3 CC1(j,k) = GG2(j,k) 1065 CC2(j,k) = GG4(j,k) 1070 continue hesshess(4,6) = zakezake(zdzd,CC1,CC2,si,nnn) do 1080 j=1,3 do 1075 k=1,3 CC1(j,k) = GG3(j,k) 1075 CC2(j,k) = GG3(j,k) 1080 continue hesshess(5,5) = zakezake(zdzd,CC1,CC2,si,nnn) do 1090 j=1,3 do 1085 k=1,3 CC1(j,k) = GG3(j,k) 1085 CC2(j,k) = GG4(j,k) 1090 continue hesshess(5,6) = zakezake(zdzd,CC1,CC2,si,nnn) do 1100 j=1,3 do 1095 k=1,3 CC1(j,k) = GG4(j,k) 1095 CC2(j,k) = GG4(j,k) 1100 continue hesshess(6,6) = zakezake(zdzd,CC1,CC2,si,nnn) hesshess(4,3) = hesshess(3,4) hesshess(5,3) = hesshess(3,5) hesshess(6,3) = hesshess(3,6) hesshess(5,4) = hesshess(4,5) hesshess(6,4) = hesshess(4,6) hesshess(6,5) = hesshess(5,6) hesshess(1,1) = hesshess(1,1) - ( (n - dfloat(msamp)) / th3) aaaa = -(yotsum - ( (n - dfloat(msamp)) * th1)) / (th3 ** 2) hesshess(1,3) = hesshess(1,3) + aaaa hesshess(3,1) = hesshess(3,1) + aaaa aaaa = (yotysum - (2 * th1 * yotsum) + 6 ( (n - dfloat(msamp)) * th1 * th1)) aaaa = -(aaaa / (th3 ** 3)) + 6 ( (n - dfloat(msamp)) / (2 * th3 * th3)) hesshess(3,3) = hesshess(3,3) + aaaa return end subroutine fullmlso(th1,th2,th3,th4,th6,th7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,n,msamp,ierrcde2,zmuxmle, 6 sigxmle,sigumle,sigemle,beta0mle,beta1mle,rhoqtmle) implicit real*8(a-h,o-z) real*8 hesshess(6,6),derder(6,1),told(6,1),tnew(6,1) real*8 tempa(6,1),tempd(6,6),tempf(6,6) epss = 100. ierrcde2 = 0 told(1,1) = th1 told(2,1) = th2 told(3,1) = th3 told(4,1) = th4 told(5,1) = th6 told(6,1) = th7 do 100 j=1,6 100 tnew(j,1) = told(j,1) c print *, " " c print *, "************************************" c print *, "fullmlse, start" c print *, "Starting values = ",(told(k,1),k=1,6) do 300 kk=1,9999 do 200 j=1,6 200 told(j,1) = tnew(j,1) t1 = told(1,1) t2 = told(2,1) t3 = told(3,1) t4 = told(4,1) t6 = told(5,1) t7 = told(6,1) c print *, "Calling fullmlhe, kk=",kk call fullmlhe(t1,t2,t3,t4,t6,t7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,hesshess,n,msamp) c print *, "Calling fullmlde, kk=",kk call fullmlde(t1,t2,t3,t4,t6,t7,m1,m3,zdtzd,zdsum, 6 yotsum,yotysum,derder,n,msamp) c print *, "Finished fullmlde, kk=",kk do 363 j=1,6 derder(j,1) = derder(j,1) / dfloat(n) do 362 k=1,6 hesshess(j,k) = -hesshess(j,k) / dfloat(n) 362 continue 363 continue hessmax = .1 do 489 j=1,6 derder(j,1) = derder(j,1) / hessmax do 488 k=1,6 488 hesshess(j,k) = hesshess(j,k) / hessmax 489 continue c print *, "Calling xinverse, kk=",kk iezcode = 0 call zinverse(hesshess,tempd,6,bigdet,iezcode) c print *, "Finished zinverse, bigdet=",bigdet if (iezcode.eq.0) bigdet = bigdet ** (1.0d0 / 6.0d0) cmin = 1.0d-50 if (dabs(bigdet).le.cmin) then call ones(told,6,1) ierrcde2 = 1 goto 400 endif if (dabs(bigdet).ge.cmin) then call zmmultm(tempd,derder,tempa,6,6,1) do 250 j=1,6 250 tnew(j,1) = told(j,1) + tempa(j,1) ggg = 0. do 251 j=1,6 251 ggg = ggg + dabs(tnew(j,1)) if (ggg.ge.100000.0d0) go to 400 c print *, " " c print *, "tnew = ",(tnew(k,1),k=1,6) endif epss = 0. do 275 j=1,6 275 epss = epss + dabs(tnew(j,1) - told(j,1)) if ((epss.le.0.0001d0).and.(kk.ge.5)) then goto 400 endif 300 continue 400 continue if ((epss.le.0.0001d0).and.(kk.ge.5)) kk = 999999 if (kk.le.9999) ierrcde2 = 1 tml1 = told(1,1) tml2 = told(2,1) tml3 = told(3,1) tml4 = told(4,1) tml6 = told(5,1) tml7 = told(6,1) rhoqtmle = tml4 / dsqrt(tml3 * tml7) zmuxmle = tml2 sigxmle = dsqrt(abs(tml7)) sigumle = dsqrt(abs(tml6 - tml7)) sigemle = dsqrt(abs(tml3 - tml4)) beta1mle = tml4 / tml7 beta0mle = tml1 - (beta1mle * tml2) if (ierrcde2.eq.1) then print *, " " print *, "************************************" print *, "MLE did not Converge" print *, "Value of ierrcde2 = ",ierrcde2 print *, "************************************" print *, " " endif return end double precision function zakezake(zdzd,AA1,AA2,si,msamp) implicit real*8(a-h,o-z) real*8 zdzd(3,3),AA1(3,3),AA2(3,3),si(3,3) real*8 bb1(3,3),bb2(3,3),bb3(3,3),bb4(3,3) c cc1 = (nnn./2) .* qtrace(si * AA1 * si * AA2); call zmmultm(si,AA1,bb1,3,3,3) call zmmultm(bb1,si,bb2,3,3,3) call zmmultm(bb2,AA2,bb3,3,3,3) sszz = trace(bb3,3) cc1 = (dfloat(msamp) / 2.0d0) * sszz c cc2 = -.5 .* qtrace(si * AA1 * si * AA2 * si * zdzd); call zmmultm(si,AA1,bb1,3,3,3) call zmmultm(bb1,si,bb2,3,3,3) call zmmultm(bb2,AA2,bb3,3,3,3) call zmmultm(bb3,si,bb2,3,3,3) call zmmultm(bb2,zdzd,bb4,3,3,3) sszz = trace(bb4,3) cc2 = -.5 * sszz c cc3 = -.5 .* qtrace(si * AA2 * si * AA1 * si * zdzd); call zmmultm(si,AA2,bb1,3,3,3) call zmmultm(bb1,si,bb2,3,3,3) call zmmultm(bb2,AA1,bb3,3,3,3) call zmmultm(bb3,si,bb2,3,3,3) call zmmultm(bb2,zdzd,bb4,3,3,3) sszz = trace(bb4,3) cc3 = -.5 * sszz zakezake = cc1 + cc2 + cc3 return end