SUBROUTINE GMM(gmm_type,a,wtmtx_a,gmm_a,weighting_type,nma,        &
               init_estimate,gmm_itts,do_stats,deriv_type,         &
               eps_vector,itmax,FRET,covmtx,use_gmm_a,indep_var,   &
               nerr)
implicit none

include  'dim_dat.inc'

!maximum possible number of moments
integer, parameter ::max_num_moments=no_of_equations*              &
                                     no_of_instruments
! the size of the first dimension of Fmtx
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)

real(8), dimension(:)                ::a,eps_vector,wtmtx_a,gmm_a
real(8), dimension(:,:)              ::covmtx
real(8), dimension(:,:)              ::indep_var
integer  weighting_type,deriv_type,nerr,gmm_itts,itmax,HESSINreset
integer  nma,init_estimate, do_stats, GMM_type
real(8)  fret

real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs)::Fmtx
real(8), dimension(max_num_moments,max_num_moments)  ::W
real(8), dimension(size(a))                          ::a1
real(8), dimension(1:1,size(a))                      ::as
real(8), dimension(no_of_equations,no_of_equations)  ::crosscovmtx
real(8), dimension(1:1,no_of_equations,nobs)         ::resid
integer, dimension(no_of_equations)                  ::num_inst
integer, dimension(3)                                ::spare_params
character(len=no_of_instruments),dimension(                        &
                                     no_of_equations)::instruments
real(8)  conv1,conv2,conv3,schlength,setconv1,setconv2,setconv3
real(8)  setgmmconv1,setgmmconv2,setgmmconv3
real(8)  small,root_small,gmm_root_small,var,gmm_itt_small
integer  i,j,ii,num_moments,num_derivs,m,gmm_counter
logical  nsls_type,gmm_stage,matrix_ill,use_gmm_a
character(len=10) nsls_string

common/flag_instruments/ instruments
common/gmm_specific/     Fmtx,W,num_inst

external funcgmm0
external funcgmm1

include 'blas.inc'
include 'routines.inc'

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

interface
   SUBROUTINE calccovmtx(a,deriv_type,eps_vector,covmtx,           &
                         num_moments,matrix_ill)
      real(8), dimension(:)             ::a,eps_vector
      real(8), dimension(:,:)           ::covmtx
      integer  deriv_type,num_moments
      logical  matrix_ill
   end SUBROUTINE calccovmtx
END INTERFACE

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

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

  if(size(eps_vector).ne.size(a)) then
    write(*,*)'incorrect dimensioning between parameter vector'
    write(*,*)'and eps_vector:  PROGRAM STOPPED'
    STOP
  end if

! force an error if certain files do not exist.
   write(*,*) 'testing the opening of file a.vec'
   open(10,file='a.vec',status='old')
   close(10)
   write(*,*) 'testing the opening of file wtmtxa.vec'
   open(10,file='wtmtxa.vec',status='old')
   close(10)
   write(*,*) 'testing the opening of file weight.in'
   open(10,file='weight.in',status='old')
   close(10)
   write(*,*) 'testing the opening of file weight.out'
   open(10,file='weight.out',status='old')
   close(10)
   write(*,*) 'testing the opening of file series.out'
   open(10,file='series.out',status='old')
   close(10)

! force an error if nma is not set to a valid option
   if(nma.lt.-3) then
      write(*,*)'invalid NMA choice - PROGRAM STOPPED'
      stop
   end if

! force an error if new GMM procedure is used with nma=-3
   if((gmm_type.eq.1).and.(nma.eq.-3)) then
      write(*,*)'incompatible NMA choice and GMM TYPE:'
      write(*,*)' PROGRAM STOPPED'
      stop
   end if

! define a small number
   small=epsilon(0.0d0)
   root_small=small**0.5d0
   gmm_root_small=small**0.5d0
   gmm_itt_small=small**0.3d0

! define parameters for minimisation routine (see MINIMIZE)
   setconv1 = root_small
   setconv2 = size(a) * root_small
   setconv3 = size(a) * root_small
   setgmmconv1 = gmm_root_small
   setgmmconv2 = size(a) * gmm_root_small
   setgmmconv3 = size(a) * gmm_root_small

   HESSINreset = 6
   schlength = 10.0d0
   nerr=0
   num_derivs=2

! set num_inst and num_moments
   num_moments=0
   num_inst=0
   do j=1,no_of_equations
     do i=1,no_of_instruments
       if (instruments(j)(i:i).eq.'1') then 
         num_moments=num_moments+1
         num_inst(j)=num_inst(j)+1
       end if
     end do
   end do

! determine if all the instruments are the same,
! if so NSLS_TYPE=.true.
   gmm_stage=.false.
   nsls_type=.true.
   nsls_string='(modified)'
   do j=2,no_of_equations
     if (instruments(j).ne.instruments(1)) nsls_type=.false.
   end do
   if(nsls_type) nsls_string=''

! parameters which must be passed to function to be minimised.
   spare_params(1)=weighting_type
   spare_params(2)=nma
   spare_params(3)=num_moments

   if(init_estimate.eq.2 .or. init_estimate.eq.3) then    ! do 2SLS
     call setinst
     m = 1
     W=0.0d0
     do j=1,no_of_equations
       i=m+num_inst(j)-1
       W(m:i,m:i) = matmul( Fmtx(1,m:i,:),transpose(Fmtx(1,m:i,:)) )
       call inverse(W(m:i,m:i),matrix_ill)
       if(matrix_ill) then
           write(*,*)'PROGRAM STOPPED'
           STOP
       end if
       m=m+num_inst(j)
     end do
     conv1 = setconv1
     conv2 = setconv2 
     conv3 = setconv3
     call MINIMIZE (spare_params,a,max_num_as,funcgmm0,            &
                    HESSINreset,itmax,schlength,deriv_type,        &
                    eps_vector,num_derivs,conv1,conv2,conv3,FRET,  &
                    nerr)
     write(*,*)'function stepsize =',conv1
     write(*,*)'L2 norm of change in parameter vector =',conv2
     write(*,*)'L2 norm of the gradient',conv3
     write(*,*)
     write(*,*)
     write(*,*) '2SLS ',nsls_string, 'parameter estimates:'
     write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
     write(*,*) 'parameter vector = ',a
     if(do_stats.eq.2) then
       call cross_equ_cov(a,crosscovmtx)
       var=0.0d0
       do i=1,no_of_equations
         var=var+crosscovmtx(i,i)
       end do
       var=var/real(no_of_equations,kind(0.0d0))
       W(1:num_moments,1:num_moments)=W(1:num_moments,             &
                                        1:num_moments)/var
       call write_stats(a,gmm_stage)
       W(1:num_moments,1:num_moments)=W(1:num_moments,1:           &
                                        num_moments)*var
     end if
   wtmtx_a=a
   end if

   if(init_estimate.eq.3) then  ! do 3SLS. Uses existing W from 2SLS
     call cross_equ_cov(a,crosscovmtx)
     if (nsls_type) then
       call inverse(crosscovmtx,matrix_ill)
       if(matrix_ill) then
          write(*,*)'PROGRAM STOPPED'
          stop
       end if     
       W(1:num_moments,1:num_moments)=mtx_cross_product(           &
                  crosscovmtx,W(1:num_inst(1),1:num_inst(1)))
     else
       m = 1
       do j=1,no_of_equations
         i=m+num_inst(j)-1
         W(m:i,m:i) = W(m:i,m:i)/crosscovmtx(j,j)
         m=m+num_inst(j)
       end do
     end if
     conv1 = setconv1
     conv2 = setconv2 
     conv3 = setconv3
     call MINIMIZE (spare_params,a,max_num_as,funcgmm0,            &
                    HESSINreset,itmax,schlength,deriv_type,        &
                    eps_vector,num_derivs,conv1,conv2,conv3,       &
                    FRET,nerr)
     write(*,*)'function stepsize =',conv1
     write(*,*)'L2 norm of change in parameter vector =',conv2
     write(*,*)'L2 norm of the gradient',conv3
     write(*,*)
     write(*,*)
     write(*,*) '3SLS ',nsls_string, 'parameter estimates:'
     write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
     write(*,*) 'parameter vector = ',a
     if(do_stats.eq.2) call write_stats(a,gmm_stage)
     wtmtx_a=a
   end if

   if(init_estimate.eq.0) then    ! do GMM with W = identity
      W(1:num_moments,1:num_moments)=IDENTITY(num_moments)
      conv1 = setconv1
      conv2 = setconv2 
      conv3 = setconv3
      call MINIMIZE (spare_params,a,max_num_as,funcgmm0,           &
                     HESSINreset,itmax,schlength,deriv_type,       &
                     eps_vector,num_derivs,conv1,conv2,conv3,      &
                     FRET,nerr)
      write(*,*)'function stepsize =',conv1
      write(*,*)'L2 norm of change in parameter vector =',conv2
      write(*,*)'L2 norm of the gradient',conv3
      write(*,*)
      write(*,*)
      WRITE(*,*) 'GMM with W=IDENTITY parameter estimates:' 
      write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
      write(*,*) 'parameter vector = ',a
      wtmtx_a=a
   end if

   if(init_estimate.eq.1) then    ! do GMM with W = specified
      open(10,file='weight.in',status='old')
      do i=1,num_moments
        do j=1,num_moments
          read(10,*)W(i,j)
        end do
      end do
      close(10)
      write(*,*)
      write(*,*)'weighting matrix used for initial estimate'
      do i=1,num_moments
         write(*,*)W(i,:)
      end do
      write(*,*)

      W(1:num_moments,1:num_moments)=W(1:num_moments,1:            &
                      num_moments)/real(nobs,kind(0.0d0))

      conv1 = setgmmconv1
      conv2 = setgmmconv2 
      conv3 = setgmmconv3
      call MINIMIZE (spare_params,a,max_num_as,funcgmm0,           &
                     HESSINreset,itmax,schlength,deriv_type,       &
                     eps_vector,num_derivs,conv1,conv2,conv3,      &
                     FRET,nerr)
      write(*,*)'function stepsize =',conv1
      write(*,*)'L2 norm of change in parameter vector =',conv2
      write(*,*)'L2 norm of the gradient',conv3
      write(*,*)
      write(*,*)
      write(*,*) 'GMM with W=specified, parameter estimates:' 
      write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
      write(*,*) 'parameter vector = ',a
      if (do_stats.eq.2) call write_stats(a,gmm_stage)
      wtmtx_a=a
   end if


   gmm_stage=.true.
   if(use_gmm_a) a=gmm_a

   if (gmm_type.eq.0) then
     a1=wtmtx_a
     gmm_counter=0
     write(*,*)'STANDARD GMM PROCEDURE (GMM TYPE 0) NOW STARTING'
     write(*,*)
     do
       gmm_counter=gmm_counter+1
       call wtmtx(wtmtx_a,weighting_type,nma,num_moments)
       conv1 = setgmmconv1
       conv2 = setgmmconv2 
       conv3 = setgmmconv3
       call MINIMIZE (spare_params,a,max_num_as,funcgmm0,          &
                      HESSINreset,itmax,schlength,deriv_type,      &
                      eps_vector,num_derivs,conv1,conv2,conv3,     &
                      FRET,nerr)
       write(*,*)'function stepsize =',conv1
       write(*,*)'L2 norm of change in parameter vector =',conv2
       write(*,*)'L2 norm of the gradient',conv3
       write(*,*)
       write(*,*)
       wtmtx_a=a
       open(10,file='wtmtxa.vec',status='old')
       rewind(10)
       do ii=1,size(wtmtx_a)
         write(10,*)wtmtx_a(ii)
       end do
       close(10)
       write(*,*) 'GMM minimisation ',gmm_counter,' complete' 
       write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
       write(*,*) 'parameter vector = ',a
       write(*,*)
       if (do_stats.eq.2) call write_stats(a,gmm_stage)
       if(norm2(a1-a).le.(gmm_itt_small*norm2(a)+small)) exit
       if(gmm_counter.eq.gmm_itts) exit
       a1=a
     end do
     if (do_stats.eq.1) call write_stats(a,gmm_stage)
   end if

   if (gmm_type.eq.1) then
     write(*,*)'NEW GMM PROCEDURE (GMM TYPE 1) NOW STARTING'
     write(*,*)
     conv1 = setgmmconv1
     conv2 = setgmmconv2 
     conv3 = setgmmconv3
     call MINIMIZE (spare_params,a,max_num_as,funcgmm1,            &
                    HESSINreset,itmax,schlength,deriv_type,        &
                    eps_vector,num_derivs,conv1,conv2,conv3,       &
                    FRET,nerr)
     write(*,*)'function stepsize =',conv1
     write(*,*)'L2 norm of change in parameter vector =',conv2
     write(*,*)'L2 norm of the gradient',conv3
     write(*,*)
     write(*,*)
     nma=spare_params(2)
     write(*,*) 'GMM complete' 
     write(*,*) 'function minimum = ',fret/real(nobs,kind(0.0d0))
     write(*,*) 'parameter vector = ',a
     write(*,*)
     if (do_stats.eq.2 .or. do_stats.eq.1) call write_stats(a,     &
                                           gmm_stage)
   end if

   open(10,file='weight.out',status='old')
   rewind(10)
   do i=1,num_moments
     do j=1,num_moments