gmm.f90

这是一个实数编码的遗传算法

! note: if M is a vector and W is a matrix, MATMUL(M,W)
! is interpreted as MTW
!         f=dot_product(matmul(M(1:num_moments),W(1:num_moments,   &
!                       1:num_moments)),M(1:num_moments))

  MTW(1:no_as,:)=matmul(M(1:no_as,:),W(1:spare_params(3),          &
                        1:spare_params(3)))
  do i=1,no_as
    F(i)=dot_product(MTW(i,:),M(i,:))
  end do

return
end SUBROUTINE funcgmm0

!  *********************************************************          
SUBROUTINE funcgmm1(spare_params,as,F,nerr)

!  Defines the GMM objective function
!  It is of the form required by MINIMIZE.

!  spare_params =  integer parameters which need to be passed.
!                  (1) = weighting_type
!                  (2) = nma
!                  (3) = num_moments
!  as           =  set of different parameter vectors; as(1,:),
!                  as(2,:) , ..., as(i,:) represents the set of
!                  parameter vectors which the
!                  objective function is calculated for.
!  F            =  Objective function value
!  nerr         =  returns a -ve value if an error occurs.

implicit none

include  'dim_dat.inc'
integer, parameter ::max_num_moments=no_of_equations*              &
                                     no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments)   ::W
integer, dimension(no_of_equations)                   ::num_inst
integer, dimension(:)                                 ::spare_params

integer  nerr,no_as,i
real(8), dimension(:,:)                            ::as
real(8), dimension(size(as,dim=1),spare_params(3)) ::M
real(8), dimension(:)                              ::F

common/gmm_specific/     Fmtx,W,num_inst

interface
  SUBROUTINE constfmtx(as)
    real(8), dimension(:,:)      ::as
  END SUBROUTINE constfmtx
END INTERFACE

interface
   SUBROUTINE wtmtx(a,weighting_type,nma,num_moments)
      real(8), dimension(:)     ::a
      integer weighting_type, nma, num_moments
   END SUBROUTINE wtmtx
END INTERFACE

!  weighting_type = spare_params(1)
!  nma            = spare_params(2)
!  num_moments    = spare_params(3)

! the sample moment vector is the average of the row vectors in FMTX
    no_as=size(as,dim=1)
    call CONSTFMTX(as)
    M(1:no_as,:)=sum(Fmtx(1:no_as,1:spare_params(3),:),dim=3)

! the GMM objective is to minimize:  transpose(M).W.M where W is
! a weighting matrix
! note: if M is a vector and W is a matrix, MATMUL(M,W) is 
!       interpreted as MTW
! f=dot_product(matmul(M(1:num_moments),W(1:num_moments,           &
!                      1:num_moments)),M(1:num_moments))

!  no need to generalise the next bit because a do
!  loop will be needed
!  when calculating W for several 'a' vectors anyway
    do i=1,no_as
      call wtmtx(as(i,:),spare_params(1),spare_params(2),          &
                         spare_params(3))
      F(i)=dot_product(matmul(M(i,:),W(1:spare_params(3),          &
                       1:spare_params(3))),M(i,:))
    end do

return
end SUBROUTINE funcgmm1

! *******************************************************

SUBROUTINE WTMTX(a,weighting_type,nma,num_moments)

! Calculate the chosen weighting matrix
!
!  a              = vector of parameters : point at which to 
!                   calculate the weighting matrix
!  weighting_type = type of weighting matrix used.
!                   If  weighting_type<=100 then the weighting
!                   matrix is based on  E[F.FT], where F is defined
!                   in Hansen-Singleton.
!                   If  weighting_type>100 then the weighting matrix
!                   is based on  E[ (F-E[F]).(F-E[F])T ], that is
!                   the expected moments are taken account of where
!                   in original Hansen-Singleton paper the E[F] is
!                   zero. Theoretically it is if there is no model
!                   mispecification! The default weighting_type = 0.
!
!                    0 or 100 = set to identity
!                    1 or 101 = Hansen & Singleton matrix
!                    2 or 102 = Newey & West (Bartlett) matrix
!                    3 or 103 = PARZEN matrix
!                    4 or 104 = QUADRATIC SPECTRAL ESTIMATOR matrix
!                    5 or 105 = TUKEY-HANNING matrix
!  nma            = number of moving-average terms in the errors;
!                   nma = -1  Newey-West [1994] Automatic Bandwith
!                             Selection
!                   nma = -2  Andrews [1991] Bandwith Selection
!                             based on modelling moment conditions
!                             generated from the initial estimates
!                             as an AR(1) process
!                   nma = -3  Andrews [1991] Bandwith Selection
!                             based on modelling moment conditions
!                             as an AR(1) process where the
!                             residuals are from the last stage in
!                             an Iterated GMM
!  num_moments     = number of moment conditions

implicit none

include  'dim_dat.inc'
integer, parameter ::max_num_moments=no_of_equations*              &
                     no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments)   ::W
integer, dimension(no_of_equations)                   ::num_inst
integer  num_moments

real(8), dimension(:)                          ::a
real(8), dimension(1:1,size(a))                ::as
real(8), dimension(num_moments,nobs)           ::mu
real(8), dimension(num_moments,num_moments)    ::Rj
real(8), dimension(num_moments)                ::rho,var
integer  weighting_type, j, nma, wttype, old_nma
real(8)  tmp,alpha1,alpha2
logical  matrix_ill,nma_flag

common/gmm_specific/     Fmtx,W,num_inst

include 'blas.inc'

interface
  SUBROUTINE CONSTFMTX(as)
    real(8), dimension(:,:)  ::as
  end SUBROUTINE CONSTFMTX
end interface

    wttype=weighting_type
    if(wttype.ge.100) wttype=wttype-100
    if((wttype.gt.5).or.(wttype.lt.1)) then
      W(1:num_moments,1:num_moments)=IDENTITY(num_moments)
      goto 10       
    end if

  old_nma=nma
  nma_flag=.false.

!!! if nma=-1 then calculate AUTOMATIC BANDWIDTH
    if(nma.eq.-1) then
      tmp=real(nobs,kind(tmp))/100.0d0
      if(wttype.eq.1) nma=int(4.0d0*(tmp**(0.25d0)))
      if(wttype.eq.2) nma=int(4.0d0*(tmp**(2.0d0/9.0d0)))
      if(wttype.eq.3) nma=int(4.0d0*(tmp**(0.16d0)))
      if(wttype.eq.4) nma=int(4.0d0*(tmp**(0.08d0)))
      if(wttype.eq.5) nma=int(4.0d0*(tmp**(0.25d0)))
    end if

    as(1,:)=a
    call constfmtx(as)

!!! if nma=-2 or -3 then calculate bandwidth using Andrews method
    if((nma.eq.-2).or.(nma.eq.-3)) then
       if(nma.eq.-3) nma_flag=.true.
       mu(:,:)=Fmtx(1,1:num_moments,:)
       do j=1,num_moments
          rho(j)=dot_product(mu(j,1:nobs-1),mu(j,2:nobs))/         &
                   dot_product(mu(j,1:nobs-1),mu(j,1:nobs-1))
          mu(j,1:nobs-1)=(mu(j,2:nobs)-rho(j)*mu(j,1:nobs-1))
          var(j)=dot_product(mu(j,1:nobs-1),mu(j,1:nobs-1))/       &
                   real(nobs-1,kind(0.0d0))
 write(*,*)'AR & VAR(innovation) for moment',j,'=',rho(j),var(j)
       end do
       mu(:,1)=4.0d0*rho*var*rho*var/((1.0d0-rho)**8)
       alpha2=sum(mu(:,1))
       mu(:,1)=4.0d0*rho*var*rho*var/(((1.0d0-rho)**6)*            &
               ((1.0d0+rho)**2))
       alpha1=sum(mu(:,1))
       mu(:,1)=var*var/((1.0d0-rho)**4)
       tmp=sum(mu(:,1))
       alpha1=alpha1*real(nobs,kind(0.0d0))/tmp
       alpha2=alpha2*real(nobs,kind(0.0d0))/tmp

!      note the -1 at the end of the following is to make
!      the Andrews method of calculating nma the same as Newey-West
!      formula which is the method adopted.
       if(wttype.eq.1) nma=nint(0.6611d0*(alpha2**(0.2d0)))-1
       if(wttype.eq.2) nma=nint(1.1447d0*(alpha1**(                &
                           0.33333333333d0)))-1
       if(wttype.eq.3) nma=nint(2.6614d0*(alpha2**(0.2d0)))-1
       if(wttype.eq.4) nma=nint(1.7462d0*(alpha2**(0.2d0)))-1
       if(wttype.eq.5) nma=nint(1.3221d0*(alpha2**(0.2d0)))-1
       write(*,*)'nma adopted =',nma
    end if

    if(weighting_type.gt.100) then ! use E[(f-Ef)(f-Ef)^T]
      mu(:,:)=spread(sum(Fmtx(1,1:num_moments,:),dim=2)/           &
                     real(nobs,kind(0.0d0)),dim=2,ncopies=nobs)
      Fmtx(1,1:num_moments,:)=Fmtx(1,1:num_moments,:)-mu(:,:)
    end if
    W(1:num_moments,1:num_moments)=matmul(Fmtx(1,1:num_moments,    &
                         :),transpose(Fmtx(1,1:num_moments,:)))

    DO j = 1, nma
      Rj=matmul(Fmtx(1,1:num_moments,j+1:nobs),transpose(          &
                Fmtx(1,1:num_moments,1:nobs-j)))
      W(1:num_moments,1:num_moments)=W(1:num_moments,1:            &
            num_moments) + weight(wttype,j,nma)*(rj+transpose(rj))
    end do

    call inverse(W(1:num_moments,1:num_moments),matrix_ill)
    if(matrix_ill) then
       write(*,*)'ill weighting matrix - PROGRAM STOPPED'
       STOP
    end if

    if(nma_flag) nma=old_nma

10  return

CONTAINS
  function weight(wttype,j,nma)
    integer wttype,j,nma
    real(8) tmp1,tmp2,tmp3,weight

    tmp1=real(j,kind(0.0d0))/real(nma+1,kind(0.0d0))
    tmp2=acos(-1.0d0)*tmp1     ! note acos(-1)=PI
    tmp3=1.2d0*tmp2

    if(wttype.eq.1) weight=1.0           ! HANSEN-SINGLETON
    if(wttype.eq.2) weight=1.0d0-tmp1    ! NEWEY-WEST (BARTLETT)
    if(wttype.eq.3) then                 ! PARZEN
      if(tmp1.le.0.5d0) weight=1.0d0-6.0d0*(tmp1**2)+6.0d0*(tmp1**3)
      if(tmp1.gt.0.5d0) weight=2.0d0*((1.0d0-tmp1)**2)
    end if
    if(wttype.eq.4) weight=((sin(tmp3)/tmp3)-cos(tmp3))*3.0d0/     &
                           (tmp3*tmp3) ! QUADRATIC SPECTRAL
    if(wttype.eq.5) weight=1.0d0+(cos(tmp2)/2.0d0) ! TUKEY-HANNING
  end function weight

end SUBROUTINE WTMTX

! *******************************************************

subroutine GETGRAD(a,grad,deriv_type,eps_vector)

! this routine calculates the gradient of the residuals
!
! note: residual must be of the form:
!                SUBROUTINE residual(as,resid,nerr)
!                     real(8), dimension(:,:)      ::as
!                     real(8), dimension(:,:,:)    ::resid
!                     integer nerr
!                end SUBROUTINE residual
!
!
!  a               = vector of parameters
!  deriv_type      = type of numerical derivative used
!                    0 = forward approximation
!                    1 = symmetric approximation
!                    2 = stabilised symmetric approximation
! eps_vector       = a vector defining "a small number"
!                    corresponding to the parameters in 'a':
!                    if the number set in eps_vector is:
!                    +ve  then it is a proportional epsilon:
!                         delta(i) = eps_vector(i)*a(i)
!                    -ve  then it is an absolute epsilon:
!                         delta(i) = -eps_vector(i)
!                     0   then the default value (machine driven)
!  grad            = stores the gradient of the residuals

implicit none
real(8), dimension(:)                   ::a,eps_vector
real(8), dimension(5*size(a),size(a))   ::as
real(8), dimension(:,:,:)               ::grad
real(8), dimension(5*size(a),size(grad,dim=2),                     &
                       size(grad,dim=3))::term
real(8), dimension(size(a))             ::eps
real(8)  small
integer  k, deriv_type,nerr,sizea

interface
  SUBROUTINE residual(as,resid,nerr)
    real(8), dimension(:,:)       ::as
    real(8), dimension(:,:,:)     ::resid
    integer nerr
  end SUBROUTINE residual
end interface

   if((deriv_type.lt.0).or.(deriv_type.gt.2)) then
     write(*,*)'error in gradient calculation'
     stop
   end if
  sizea=size(a)
  small = sqrt(epsilon(0.0d0))

! set up as
   as=spread(a,dim=1,ncopies=sizea*5)
   where(eps_vector>0.0d0) eps=max(eps_vector*abs(a),small)
   where(eps_vector<0.0d0) eps=max(-eps_vector,small)
   where(eps_vector.eq.0.0d0) eps=max(small*abs(a),small)

   DO k=1,sizea
     as(sizea+k,k)=a(k)-eps(k)
     as(2*sizea+k,k)=a(k)+eps(k)
     as(3*sizea+k,k)=a(k)-2.0d0*eps(k)
     as(4*sizea+k,k)=a(k)+2.0d0*eps(k)
   end do

   if(deriv_type.eq.0) then
     CALL residual(as(1:3*sizea,:),term(1:3*sizea,:,:),nerr)
     DO k=1,sizea
       GRAD(k,:,:) = (term(2*sizea+k,:,:)-term(k,:,:))/EPS(k)
     end do
   end if

   if(deriv_type.eq.1) then
     CALL residual(as(sizea+1:3*sizea,:),term(sizea+1:3*sizea,     &
                                                   :,:),nerr)
     DO k=1,sizea
       GRAD(k,:,:)=(term(2*sizea+k,:,:)-term(sizea+k,:,:))/        &
                                            (2.0d0*EPS(k))
     end do
   end if

   if(deriv_type.eq.2) then
     CALL residual(as(sizea+1:5*sizea,:),term(sizea+1:5*sizea,     &
                                                   :,:),nerr)
     DO k=1,sizea
       GRAD(k,:,:)=(term(3*sizea+k,:,:)-term(4*sizea+k,:,:) +      &
              8.0d0*(term(2*sizea+k,:,:)-term(sizea+k,:,:)))/      &
              (12.0d0*EPS(k))
     end do
   end if

   if(nerr.lt.0) then
     write(*,*)'error occurred in RESIDUAL calculation'
     write(*,*)'PROGRAM STOPPED'
     stop
   end if

RETURN
END subroutine GETGRAD