c
c     Computes portmanteau statistic for multivariate time series
c     regression for a ik-dimensional process.
c
c    The program has been checked, but some errors may exist. 
c    It is intended for use in my courses only. Ruey S. Tsay, December 1990.
c
      implicit double precision(a-h,o-z)
      parameter(mxn=100000, id=10, mlag = 120)

      dimension e(mxn,id), cov(id,id), corr(id,id), 
     *     trcorr(id,id), sigma(id,id), esq(mxn,id)
      dimension siginv(id,id) 

c     Read in data
c
      call setup(e,mxn,nob,id,ik,nlag,npar)
c
c  remove sample means
c
      do i=1, ik
       ave  = 0.0d0
       do  j=1, nob
        ave = ave + e(j,i)
       enddo
       ave = ave/dfloat(nob)
       do j=1, nob
        e(j,i) = e(j,i) - ave
       enddo
      enddo
c
      write(6,41)
 41   format(1x,'Compute Q-stat for the squared series? (1=y): ',$)
      read(5,*) isq
c
      if(isq .ne. 1)go to 45
c------------- square series
      do i=1, ik
       do it=1, nob
        esq(it,i) = e(it,i)*e(it,i)
       enddo
      enddo
c
c  remove sample means
c
      do i=1, ik
       ave  = 0.0d0
       do  j=1, nob
        ave = ave + esq(j,i)
       enddo
       ave = ave/dfloat(nob)
       do j=1, nob
        esq(j,i) = esq(j,i) - ave
       enddo
      enddo
c
 45   continue
      print*,'Output file: fort.36'

      tt = dlog(dfloat(nob))
      idlag=tt

      write(36,*)'Suggested # of Lags = ',idlag
c      print*,'Idea # of lags = ',idlag
c
c     Compute sample covariance matrix

      call covar(e,sigma,mxn,nob,id,ik,0)
c
      call mtinv(sigma,siginv,id,ik)
c
      q = 0.0d0
      qcross = 0.0d0

      iend = nlag
      if(iend .gt. mlag)iend= mlag
c
      write(36,38)
 38   format(2x,'m',7x,'Orig-Q',7x,'LM-Q',7x,'LB-Q',5x,'d.f.')
      do it = 1, iend
        call covar(e,cov,mxn,nob,id,ik,it)
        call correlation(cov,corr,trcorr,siginv,id,ik)
        call mult(corr,trcorr,trace,trcross,mxn,nob,id,ik,it)
        q = q + trace
        qcross = qcross + trcross
        qstar = q + dfloat((ik**2)*it*(it+1))/(2.0d0*dfloat(nob))
        idf = it*ik*ik-npar
        if(idf .lt. 0)idf = 0
c
        write(36,37) it, q, qstar, qcross, idf
      enddo
   37   format(1x,i3,1x,3f12.4,i5)
cc****
c*** squared series
c****
      if(isq .ne. 1)stop
       write(36,42)
 42    format(1x,'*** For squared series *** ')
c
c     Compute sample covariance matrix of the squared series

      call covar(esq,sigma,mxn,nob,id,ik,0)
c
      call mtinv(sigma,siginv,id,ik)
c
      q = 0.0d0
      qcross = 0.0d0
c
      write(36,38)
      do it = 1, iend
        call covar(esq,cov,mxn,nob,id,ik,it)
        call correlation(cov,corr,trcorr,siginv,id,ik)
        call mult(corr,trcorr,trace,trcross,mxn,nob,id,ik,it)
        q = q + trace
        qcross = qcross + trcross
        qstar = q + dfloat((ik**2)*it*(it+1))/(2.0d0*dfloat(nob))
        idf = it*ik*ik-npar
        if(idf .le. 0)idf = 0
c
        write(36,37) it, q, qstar, qcross, idf
      enddo
c
      stop
      end

c------------------------------------------------------------------
      subroutine covar(resid,c,n,nob,id,ik,lag)
      implicit double precision(a-h,o-z)
c
c     Computes covariance matrices

      dimension resid(n,id), c(id,id)

      do i = 1,ik
        do j = 1,ik
          c(i,j) = 0.0d0
        enddo
      enddo

      do i = 1,ik
        do j = 1,ik
          accu=0.0d0
          do k = lag + 1, nob
            accu = accu + resid(k,i)*resid(k-lag,j)
          enddo
          c(i,j) = accu
        enddo
      enddo      

c     Divide matrix by number of observations, n

      do i = 1,ik
        do j = 1,ik
          c(i,j) = c(i,j)/dfloat(nob)
        enddo
      enddo

      return
      end

c------------------------------------------------------------------

      subroutine correlation(c,r,rt,sinv,id,ik)
      implicit double precision(a-h,o-z)
c
c     Computes correlation matrices

      dimension c(id,id), r(id,id), rt(id,id), sinv(id,id)

c     Compute r = Rk

      do i = 1,ik
        do j = 1,ik
          r(i,j) = 0.0d0
        enddo
      enddo

      do i = 1,ik
        do j = 1,ik
          do k = 1,ik
c            r(i,j) = r(i,j) + sinv(i,k)*c(k,j)
            r(i,j) = r(i,j) + c(i,k)*sinv(k,j)
          enddo
        enddo
      enddo

c     Compute rt = R-k

      do i = 1,ik
        do j = 1,ik
          rt(i,j) = 0.0d0
        enddo
      enddo

      do i = 1,ik
        do j = 1,ik
          do k = 1,ik 
c            rt(i,j) = rt(i,j) + sinv(i,k)*c(j,k)
            rt(i,j) = rt(i,j) + c(k,i)*sinv(k,j)
          enddo
        enddo
      enddo

      return
      end

c------------------------------------------------------------------

      subroutine mult(r,rt,sum,sumcr,n,nob,id,ik,lag)
c
      implicit double precision(a-h,o-z)
c
c     Multiplies r and rt, computes trace, then multiplies by n.
c
      parameter(m=10)
      dimension r(id,id), rt(id,id), temp(m,m)

c     Multiply r and rt:

      do i = 1,ik
        do j = 1,ik
          temp(i,j) = 0.0d0
        enddo
      enddo

      do i = 1,ik
        do j = 1,ik
          do k = 1,ik
c            temp(i,j) = temp(i,j) + r(i,k)*rt(k,j)
            temp(i,j) = temp(i,j) + rt(i,k)*r(k,j)
          enddo
        enddo
      enddo

c     Compute trace

      sum = 0.0d0

      do i = 1,ik
        sum = sum + temp(i,i)
      enddo

      sumcr = (dfloat(nob*nob))*(sum/(dfloat(nob-lag)))

      sum = sum*dfloat(nob)

      return
      end

c---------------------------------------------------------------

      subroutine setup(resid,n,nob,id,k,nlag,npar)
      implicit double precision(a-h,o-z)

c     Read in residuals

      dimension resid(n,id), jdx(10)
      character*40  fname
c
      write(6,1)
    1 format(1x,'residual (data) file_name: ',$)
      read(5,2) fname
    2 format(a40)

      write(6,3)
    3 format(1x,'#(columns) in the data file: ',$)
      read(5,*) nc
      write(6,34)
 34   format(1x,'#(parameters) if a model was fit: ',$)
      read(5,*) npar
      write(6,35)
 35   format(1x,'maximum number of CCM mtx: ',$)
      read(5,*) nlag
      if(nlag.le.0)nalg = 12
      if(npar.le.0)npar = 0
c
      open(21,file=fname,status='old')

      nob = 0
      do 10 i = 1,n
        read (21,*,end=11) (resid(i,j), j=1,nc)
         nob = nob+1
 10     continue
 11     close(21)
c
       print*,'nob = ', nob
c
      write(6,4)
 4    format(1x,'The number of series: ',$)
      read(5,*) k

      if(k.eq.nc)return
      write(6,5) k
 5    format(1x,'Input',i3,' locators of residual series: ',$)
      read(5,*) (jdx(i),i=1,k)
c
      do 15 i=1, nob
       do 13 j=1, k
 13       resid(i,j) = resid(i,jdx(j))
 15       continue
c
      return
      end
C*******************************************************************
      SUBROUTINE MTINV(A,DA,IDM,IDIM)
C***********************************************************************
C     MTINV COMPUTES THE INVERSE MATIRX OF A
C**********************************************************************
C**** 
      implicit double precision(a-h,o-z)
      dimension A(IDM,IDM),DA(IDM,IDM)
C****
      DETA=1.0D0
      IF(IDIM .EQ. 1) GO TO 600                                    
C****                                                                         
      DO 100 I=1,IDIM                                                         
      DO 100 J=1,IDIM                                                         
  100 DA(I,J)=A(I,J)                                                          
C****                                                                         
  120 DO 500 I=1,IDIM                                                         
      PIVOT=DA(I,I)                                                           
      DETA=DETA*PIVOT                                                         
C****                                                                         
C**** DIVIDE PIVOT ROW BY PIVOT ELEMENT.                                      
C****                                                                         
      DA(I,I)=1.0D0
      DPIVOT=PIVOT+1.0D-20                                              
      DPIVOT=DA(I,I)/DPIVOT                                             
      PIVOT=DPIVOT                                                      
      DO 200 J=1,IDIM                                                         
  200 DA(I,J)=DA(I,J)*PIVOT                                                   
C****                                                                         
C**** REDUCE NON-PIVOT ROWS                                                   
C****                                                                         
  210 DO 500 II=1,IDIM                                                        
      IF(II.EQ.I) GO TO 500                                                   
      T=DA(II,I)                                                              
      DA(II,I)=0.0D0
      DO 300 J=1,IDIM                                                         
  300 DA(II,J)=DA(II,J)-DA(I,J)*T                                             
C****                                                                         
  500 CONTINUE                                                          
      RETURN  
  600 DA(1,1)=DETA/A(1,1)
      RETURN
      END                                                                     
C****                                                                         
C**** END OF 'MTINV'                                                          
C****                                                                         
C***********************************************************************      
