我写了一个程序，主程序是 test.f90, 内容为
!include "sim.f90"

program test
!use sim_int
implicit none
integer, parameter :: n=1000
integer, parameter :: p=2,pgarch=1,qgarch=1
real*8 :: tol=1.0d-10
real*8 :: ar_para(2,p+1),garch_para(2,1+pgarch+qgarch), beta(p),paraft(2,2)
real*8 :: ar_para1(2,p+1),garch_para1(2,1+pgarch+qgarch), beta1(p),paraft1(2,2)
real*8 :: y(n),h(n),eps(n),eta(n),eta2(n), pstay(n),tp(n,2,2)
real*8 :: cp(n,2)
real*8 :: sjsp(n,2,2),ssp(n,2)
real*8 :: eta_s(n,2),h_s(n,2)
real*8 :: x(n,1+p)
real*8 :: sp(2)
integer :: rgm(n),i,t

real*8 :: dhsg(n,2,2,1+qgarch+pgarch)
real*8 :: dhg(n,2,1+qgarch+pgarch)
real*8 :: dlsg(n,2,2,1+qgarch+pgarch)
real*8 :: dlg(n,2,1+qgarch+pgarch)

real*8 :: desa(n,2,2,1+p)
real*8 :: dea(n,2,1+p)
real*8 :: dhsa(n,2,2,1+p)
real*8 :: dha(n,2,1+p)
real*8 :: dlsa(n,2,2,1+p)
real*8 :: dla(n,2,1+p)
real*8 :: of

real*8, external :: ft, fp11
external sim,cal_cp, cal_eta_h, cal_ssp, cal_tp, cdap, cdgp  
external cof
external ugp
!ar_para(1,:)=(/0.0d0, 0.3d0,0.5d0/)
ar_para(1,:)=0.0d0
!ar_para(2,:)=(/0.0d0, 1.0d0, 0.0d0/)
ar_para(2,:)=0.0d0
garch_para(1,:)=(/0.4d0, .4d0,.5d0/)
garch_para(2,:)=(/0.1d0, .2d0,.2d0/)
beta=1.0d0
paraft(1,:)=(/1.0d0,2.0d0/)
paraft(2,:)=(/-1.0d0,2.0d0/)


call sim( n, p, qgarch, pgarch, ar_para, garch_para, ft, beta, paraft,&
         y, rgm, h, eps, eta,pstay)

x(:,1)=1.0d0
do i=1,p
        x(1:i,i+1)=sum(y)/n
        x(i+1:n,i+1)=y(1:n-i)
end do
!do i=n-50,n
!print *,x(i,:)
!end do
!pause
call cal_tp(n,p,y,ft,beta,paraft,tp)

open(1,file="tp.dat",status="replace")
write(1,fmt='(4F15.8)') (tp(t,:,:),t=1,n)
close(1)
!initial ssp estimates
ssp(1,:)=0.5d0
do t=2,n
        ssp(t,:)=matmul(ssp(t-1,:),tp(t,:,:))
        !print *, ssp(t,:)
        !if(t.eq.100) pause
end do

call cal_eta_h(n,p,pgarch,qgarch, x, y, &
                ar_para, garch_para,ssp, &
                eta_s,eta,eta2,h_s,h)

do i=1,2
        call cal_cp(n,eta_s(:,1),sqrt(h_s(:,1)),cp(:,i))
end do

sp=0.5d0
call cal_ssp(n,cp,tp,sp,sjsp,ssp)
!do t=1,n
!write(*, '(I2,5F15.8)') rgm(t), pstay(t),tp(t,1,1),tp(t,2,2), ssp(t,1), ssp(t,2
)
!end do

call cdgp(n,qgarch,pgarch,&
         ssp,eta_s,eta2,h_s,h,garch_para,&
         dhsg,dhg,dlsg,dlg)

!do t=1,30
!write(*, '(I4, 15(D15.6))') t, (dlg(t,i,:),i=1,2)
!end do
call cdap(n,p,qgarch,pgarch,&
        x,ssp,eta_s,eta,eta2,h_s,h,garch_para,&
        desa,dea,dhsa,dha,dlsa,dla)
 
call cof(n,ssp,eta_s,h_s,of)
print *,"of=",of
call ugp(n,p,pgarch,qgarch, x, y, &
        ar_para, garch_para,ssp,dlg,tol, &
        garch_para1,of)
print *,garch_para1
write(*,'(A)') "Para for GARCH:"
write(*,'(A,20D20.13)') "Para for R1=",garch_para(1,:)
write(*,'(A,20D20.13)') "Para for R2=",garch_para(2,:)
write(*,'(A)') " estimated value:"
write(*,'(A,20D20.13)') "Para for R1=",garch_para1(1,:)
write(*,'(A,20D20.13)') "Para for R2=",garch_para1(2,:)

call ur1tpp(n,p,y,sjsp(:,1,1),beta,paraft(1,:),&
        fp11,ft,tol,&
        paraft1(1,:))
write(*,'(A,2D20.13)') " paraft in R1=",paraft(1,:)
write(*,'(A,2D20.13)') "estimated val=",paraft1(1,:)

end program test

其中需要调用 ugp， 我把它写在另外一个文件 ugp.f90 中：
        !ugp= Update Garch Parameters
        ! This subroutine update the GARCH parameters in the MS_GARCH model.
        subroutine ugp(n,p,pgarch,qgarch, x, y, &
                ar_para, garch_para,ssp,tol, &
                garch_para1,of)
        use linear_operators
        implicit none
        integer,intent(in) :: n
        integer,intent(in) :: p
        integer,intent(in) :: qgarch
        integer,intent(in) :: pgarch
        real*8,intent(in) :: x(n,p+1)
        real*8,intent(in) :: y(n)
        real*8,intent(in) :: ar_para(2,1+p)
        real*8,intent(in) :: garch_para(2,1+qgarch+pgarch)
        real*8,intent(in) :: ssp(n,2)
        !real*8,intent(in) :: dlg(n,2,1+qgarch+pgarch)
        real*8,intent(in) :: tol

        real*8,intent(out) :: garch_para1(2,1+qgarch+pgarch)
        real*8,intent(out) :: of

        !local variables
        real*8 :: dlg(n,2,1+qgarch+pgarch)
        real*8 :: m1(n,2*(1+qgarch+pgarch))
        real*8 :: b(2*(1+qgarch+pgarch))
        real*8 :: v(2*(1+qgarch+pgarch))
        real*8 :: a(2*(1+qgarch+pgarch),2*(1+qgarch+pgarch))
        real*8 :: ai(2*(1+qgarch+pgarch),2*(1+qgarch+pgarch))
        integer :: i,t,errorflag,k
        real*8 :: of0

        real*8 :: eta_s(n,2)
        real*8 :: eta(n)
        real*8 :: eta2(n)
        real*8 :: h_s(n,2)
        real*8 :: h(n)
        real*8 :: g_para0(2*(1+qgarch+pgarch))
        real*8 :: g_para1(2*(1+qgarch+pgarch))
        real*8 :: dif
        external cof
        external cal_eta_h
        external cdgp
        !external findinv
      
        garch_para1=garch_para 
        !calculate initial value of of0
        call cal_eta_h(n,p,pgarch,qgarch, x, y, &
                ar_para, garch_para1,ssp, &
                eta_s,eta,eta2,h_s,h)
        call cof(n,ssp,eta_s,h_s,of0)
        print *,"of0=",of0

        do k=1,1000
                !write(*,'(A,I3)') "ITERAION=",k
                !calculate DLG
                call cal_eta_h(n,p,pgarch,qgarch, x, y, &
                        ar_para, garch_para1,ssp, &
                        eta_s,eta,eta2,h_s,h)
                call cdgp(n,qgarch,pgarch,&
                        ssp,eta_s,eta2,h_s,h,garch_para1,&
                        dlg)
                !find direction
                m1(:,1:1+qgarch+pgarch)=dlg(:,1,:)
                m1(:,2+qgarch+pgarch:2*(1+qgarch+pgarch))=dlg(:,2,:)
                b=sum(m1,dim=1)
                a=matmul(transpose(m1),m1)
 
 
                v=-a.ix.b
                !print *, "v using imsl="
                !write( *,'("v=",//,30F10.4)') v
                garch_para1(1,:)=garch_para1(1,:)+v(1:1+qgarch+pgarch)
                garch_para1(2,:)=garch_para1(2,:)+v(2+qgarch+pgarch:)
                call cal_eta_h(n,p,pgarch,qgarch, x, y, &
                        ar_para, garch_para1,ssp, &
                        eta_s,eta,eta2,h_s,h)
                call cof(n,ssp,eta_s,h_s,of)
                dif=of-of0
                !write(*,'(A,D20.5)') "Of=",of

                !halve steps to find next estimates
                do i=1,100
                        !print *, "i=",i
                        if(dif .le. 0.0d0) exit
                        v=v/2.0d0
                        garch_para1(1,:)=garch_para1(1,:)-v(1:1+qgarch+pgarch)
                        garch_para1(2,:)=garch_para1(2,:)-v(2+qgarch+pgarch:)

                        call cal_eta_h(n,p,pgarch,qgarch, x, y, &
                                ar_para, garch_para1,ssp, &
                                eta_s,eta,eta2,h_s,h)
                        call cof(n,ssp,eta_s,h_s,of)
                        dif=of-of0
                        !write(*,'(A,D15.6)') "DIF=", dif
                        !write(*,'(A,D20.5)') "Of=",of
                end do
                !dif=of-of0
                write(*,'(A,I3,A,D15.6,A,D15.6)') "ITERAION=",k," OF=",of, " DIF
=", dif
                write(*,'(A,20D20.13)') "Para for R1=",garch_para1(1,:)
                write(*,'(A,20D20.13)') "Para for R2=",garch_para1(2,:)

                if(dif.gt. -tol) then
                        if(dif.ge.0.0d0) then
                                print *,"CANNOT DECREASE O.F. ANY MORE, EXIT."
                        else
                                print *,"CONVERGENCE ACHEIVED, EXIT."
                        end if

                        exit
                else
                        of0=of
                        !garch_para1=garch_para1
                end if
        end do
        write(*,'(A)') "FINAL ESTIMATED RESULT OF GARCH PARAMETER:"
        write(*,'(A,20D20.13)') "Para for R1=",garch_para1(1,:)
        write(*,'(A,20D20.13)') "Para for R2=",garch_para1(2,:)
        return

        end subroutine ugp

程序结果为
……
ITERAION= 39 OF=   0.475281D+03 DIF=  -0.170530D-12
Para for R1= 0.6736895762196D-01 0.3988482251945D+00 0.6037354711764D+00
Para for R2= 0.5235228325052D-01 0.4987621074196D+00 0.4655789554840D+00
ITERAION= 40 OF=   0.475281D+03 DIF=   0.000000D+00
Para for R1= 0.6736895762196D-01 0.3988482251945D+00 0.6037354711764D+00
Para for R2= 0.5235228325052D-01 0.4987621074196D+00 0.4655789554840D+00
ITERAION= 41 OF=   0.475281D+03 DIF=   0.113687D-12
Para for R1= 0.6736895762196D-01 0.3988482251945D+00 0.6037354711764D+00
Para for R2= 0.5235228325052D-01 0.4987621074196D+00 0.4655789554840D+00
 CANNOT DECREASE O.F. ANY MORE, EXIT.
FINAL ESTIMATED RESULT OF GARCH PARAMETER:
Para for R1= 0.6736895762196D-01 0.3988482251945D+00 0.6037354711764D+00
Para for R2= 0.5235228325052D-01 0.4987621074196D+00 0.4655789554840D+00
   475.280502374522       0.000000000000000E+000  0.000000000000000E+000
  0.000000000000000E+000  0.000000000000000E+000  0.000000000000000E+000
Para for GARCH:
Para for R1= 0.4000000000000D+00 0.4000000000000D+00 0.5000000000000D+00
Para for R2= 0.1000000000000D+00 0.2000000000000D+00 0.2000000000000D+00
 estimated value:
可见传到主程序中的garch_para1 值发生了变化，请问这是为什么呢？一般的原因有哪些？