program testvar
      parameter(mt=4,ndy=9,ny=2007)
c************************************
c  Generates VAR estimates for a class of models
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+ cim2*xm(t-2)
c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Index i for the ith equation with m*p variates
c====================================
c
c  Parameters in control file Var.cnl. Use spaces between entries!
c  Place comments on the lines of cnl file after the symbol *
c
c===================================================================
c
c  Sample size=!n
c
c   n must be larger than 2*m*m*p
c
c===========================================
c
c  No. of equations=!m
c
c   Equation in VAR(p) where !m <= 20
c
c===========================================
c
c  No. of lags=!p
c
c   Lags in VAR(p)
c
c===========================================
c
c  Forecast=!q
c
c   Forecast !q steps ahead using the estimated parameters of VAR(p)
c
c===========================================
c
c  Standard deviation minimum=!mxsig
c  Standard deviation maximum=!mnsig
c
c   Maximum & minimum standard deviation for the equation errors
c   Each equation has an additive error time series (innovations). The
c   standard deviation of these errors are generated by a random draw
c   from a pseudo-random uniform (mnsig,mxsig).
c
c===========================================
c
c  Maximum damping=!rho
c  Minimum period==!per
c
c  Maximum damping value & minimum period for the system eigenvalues
c
c===========================================
c
c  Noise model=
c   After the = sign type one of the five model types:
c   gaussian, exp (exponential), double exp (double tailed exponential),
c   uniform, log normal
c
c===========================================
c
c  Switches for forecasting and getting the impulse response
c   Type 'true' on a line if you want to the use the true coefficients
c   for forcasting. If there is no character string 'true' on a line
c   then the estimated coefficients will be used in the forecasting.
c   Type 'impulse' if you want to output the impulse resonses of each
c   equation in an file called Impulses.out. If there is no character
c   string 'impulse' on a line then this file will not be opened.
c
c===========================================
c
c  The var.cnl file parameter reads above is followed by the line
c  # * end of var.cnl reads
c
c============================
      allocatable ::a,c,lam,v,z
      real a(:),mxsig,mnsig,rho,per
      real*8 c(:)
      complex lam(:),v(:),z(:)
      integer p,q
      character*700 name
      character*80 buf
      character*50 charead
      character*25 model
      character*20 dt,par
      character*2 delim
      character ch,imp,tr
      logical exs
c************************************
      inquire(file='var.cnl',exist=exs)
       if(.not.exs) then
        write(6,*) 'Control file var.cnl does not exist'
        stop
       endif
      open(1,file='var.cnl',err=1,iostat=io1,status='old')
      open(3,file='Var.out',err=2,iostat=io2)
c  Start time
      call cpu_time(gs)
c--------
      ip=0
      kt=0
      ia=0
      dowhile(kt<=17.and.ia==0)
       kt=kt+1
c  Clear buf
       do n=1,80
        buf(n:n)=' '
       enddo
c  Read line
       read(1,'(a80)',end=3,err=4,iostat=io3) buf
       ia=index(buf,'#')
       ig=index(buf,'*')
       if(ig==0) then
        write(3,*)'Place comment marker * after each line in Var.cnl'
        write(6,*)'Place comment marker * after each line in Var.cnl'
        write(3,'(a70)')buf
        write(6,'(a70)')buf
        stop
       endif
       ig=ig-1
       name(ip+1:ip+ig)=buf(:ig)
       ip=ip+ig
      enddo
      delim(1:2)='='
c  Sample size, no. of equations, no. of lags, forecast
      par='ize'
      n=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,'(/'' Error in read of sample size '',i5)')n
       write(6,'(/'' Error in read of sample size '',i5)')n
       stop
      endif
      par='tions'
      m=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,'(/'' Error in read of no. of equation '',i3)')m
       write(6,'(/'' Error in read of no. of equation '',i3)')m
       stop
      endif
      par='ags'
      p=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,'(/'' Error in read of no. of lags '',i5)')p
       write(6,'(/'' Error in read of no. of lags '',i5)')p
       stop
      endif
      par='ecast'
      q=numread(par,name,delim,ia,ier,io)
      if(ier>0) then
       write(3,'(/'' Error in read of forecast steps '',i5)')q
       write(6,'(/'' Error in read of forecast steps '',i5)')q
       stop
      endif
c  Traps
      if(m<1) then
       write(3,*)'No. of equations ',m,' < 1'
       stop
      endif
      if(m>20) then
       write(3,*)'No. of equations ',m,' > 20'
       stop
      endif
      if(2*m*m*p>=n) then
       write(3,*)'2*m*m*p ',2*m*m*p,' >= sample size used ',n
       stop
      endif
      if(p<1) then
       write(3,*)'No. of lags ',p,' < 1'
       stop
      endif
      if(q<1) then
       write(3,*)'No. of steps ahead ',q,' < 1'
       stop
      endif
      if(q>50) then
       write(3,*)'No. of steps ahead ',q,' > 50'
       stop
      endif
c  Maximum & minimum sigs for equation errors
      par='inimum'
      mnsig=rnumread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,'(/'' Error in read of minimum sig '',f7.3)')mnsig
       write(6,'(/'' Error in read of minimum sig '',f7.3)')mnsig
       stop
      endif
      par='aximum'
      mxsig=rnumread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,'(/'' Error in read of maximum sig '',f7.3)')mxsig
       write(6,'(/'' Error in read of maximum sig '',f7.3)')mxsig
       stop
      endif
      if(mxsig<=mnsig) then
       write(io,*)'Max sig for errors ',mxsig,' < min sig ',mnsig
       stop
      endif
      if(mxsig<1.e-17) then
       write(io,*)'Max sig for errors ',mxsig,' < 0'
       stop
      endif
c  Max damping & max period for system eigenvalues
      par='mping'
      rho=rnumread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,'(/'' Error in read of maximum damping '',f9.3)')rho
       write(6,'(/'' Error in read of maximum damping '',f9.3)')rho
       stop
      endif
      par='eriod'
      per=rnumread(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,'(/'' Error in read of maximum period '',f9.3)')per
       write(6,'(/'' Error in read of maximum period '',f9.3)')per
       stop
      endif
c  Trap values
      if(rho>=1.) then
       write(io,*)'Max damping term ',rho,' > 1'
       stop
      endif
      if(per<=0.0) then
       write(io,*) 'Max period ',per,' <= 0'
       stop
      endif
c-------
c  Model type
      par='model'
      dt=charead(par,name,delim,ia,ier,3)
      if(ier>0) then
       write(3,*)'Error in reading model type'
       write(6,*)'Error in reading model type'
       stop
      endif
      if(ia>0) then
        ig=index(name,'gauss')
        id=index(name,'aile')
        ie=index(name,'exp')
        iu=index(name,'unif')
        iln=index(name,'log')
        if(ig>0) then
         ch='g'
         model='gaussian'
          elseif(id>0) then
           ch='d'
            model='double exp'
             elseif(ie>0) then
              ch='e'
              model='exponential'
               elseif(iu>0) then
              ch='u'
              model='uniform'
               elseif(iln>0) then
                ch='l'
                model='log normal'
                 else
                  write(io,*) 'Model control letter is incorrect ',dt
                  stop
        endif
         elseif(ia==0) then
          write(io,'(/7x,''Error in read of innovations model type'')')
          write(6,'(/7x,''Error in read of innovations model type'')')
         stop
      endif
c  Flag for using true coefficients in forecast
      imp='n'
      tr='f'
      ig=index(name,'true')
      if(ig>0) then
       tr='t'
      endif
c  Flag for impulse response
      ig=index(name,'imp')
      if(ig>0) then
       open(7,file='Impulses.out',err=5,iostat=io5)
       imp='i'
      endif
c------
      mn=m*n
      mp=m*p
      mmp=mp*m
      ma=mp-1
c------
      i2=mn
      i3=i2+mmp
      if(m<=40) then
       i4=i3+m*n
        else
         i4=i3+1
      endif
      i5=i4+mmp
      if(m<=40) then
       i6=i5+m
        else
         i6=i5+1
      endif
      i7=i6+m
      i8=i7+m*q
      ia=i8+m*q
c------
c  Allocate space
      allocate(a(0:ia),c(0:mmp-1),lam(0:ma),z(0:ma),v(0:4*mp+14),stat=i
     *bad)
      if(ibad/=0) then
       ntot=ia+2*(mmp+2*ma+4*mp+15)
       write(3,*)'Unable to allocate work space for arrays a & c'
       write(3,*)'Bytes needed ',ntot
       stop
      endif
c----------
      call run(a(0),a(i2),a(i3),a(i4),a(i5),a(i6),a(i7),a(i8),c,lam,v,z,
     *mnsig,mxsig,rho,per,ch,imp,tr,model,m,n,p,q,mp,3)
c+++++++++++++++++++++++++++++++++++++++++++
c  Time and date of run & run time
      call head(mt,ndy,ny,gs,rt,3)
c+++++++++++++++++++++++++++++++++++++++++++
      deallocate(a,c,lam,v,z)
      stop
c---------
  1   write(6,*) '**** Error on opening control file 1 ',io1
      stop
  2   write(6,*) '**** Error on opening output file ',io2
      stop
  3   write(3,*)'End of file on var.cnl reads'
      write(6,*)'End of file on var.cnl reads'
  4   write(3,*)'Place the symbol # after the var.cnl reads ',io3
      write(6,*)'Place the symbol # after the var.cnl reads ',io3
      stop
  5   write(6,*) '**** Error on opening impulse response file ',io5
      stop
      end
c
c===========================================
c
      subroutine run(x,b,r,s,r2,sig,e,y,c,lam,v,z,mnsig,mxsig,rho,per,ch
     *,imp,tr,model,m,n,p,q,mp,io)
c********************************************************
c  Open data file and read data
c========================================================
      integer p,q,t,seed
      real x(m*n),b(m,mp),r(m*n),s(m*mp),r2(m),sig(m),e(m*q),y(m*q),rho,
     *per,mxsig,mnsig
      real*8 c(m,mp)
      complex lam(0:mp-1),z(0:mp-1),v(4*mp+15)
      character*25 buf,model
      character ch,init,imp,tr
      intent(in) mnsig,mxsig,rho,per,ch,imp,tr,model,m,n,p,q,mp,io
c********************************************
c  Initialize generator
      seed=1234567
      call random_seed(seed)
      call random_number(r(1:m))
c  Standard deviations for equation errors
      range=mxsig-mnsig
      do k=1,m
       sig(k)=mnsig+r(k)*range
      enddo
c------
      write(io,'(''Sample size = '',i9,2x,''No.of equations = '',i3,2x,'
     *'Lags = '',i2)') n,m,p
      write(io,'(7x,''Model Type is: '',a20)') model
      write(io,'(7x,''Maximum Damping for System Eigenvalues = '',f7.4)'
     *)rho
      write(io,'(7x,''Minimum Period for System Eigenvalues = '',f9.2)')
     *per
      write(io,'(7x,''Maximum Standard Deviation for Model Errors = '',f
     *7.2)') mxsig
      write(io,'(7x,''Minimum Standard Deviation for Model Errors = '',f
     *7.2)') mnsig
      if(tr=='t') then
       write(io,'(/2x,''True VAR coeffients used in the forecasts''/)')
      endif
c----------
      n2=mp/2
      nran=n2+1
      call random_number(r(1:nran))
c  Frequency of eigenvalue
      do k=1,nran
       r(k)=r(k)/per
      enddo
      write(io,'(50(''=''))')
      write(io,'(11x,''Complex System Eigenvalues'')')
      write(io,'(3x,''Amplitudes'',7x,''Periods''/)')
      do k=1,n2
       s(k)=((-1)**(k-1))*(1.-float(k-1)/float(n2))*rho
       write(io,'(1x,i2,2x,f8.5,6x,f12.2)') k,s(k),1./r(k)
      enddo
      write(io,'(50(''=''))')
c----------
c  Eigenvalues of system
      call sffti(mp,v)
      pi=acos(-1.0)
      lam(0)=cmplx(rho,0.)
      if(2*n2==mp) then
c  mp is even
       lam(n2)=cmplx(rho,0.)
c  Rest of eigenvalues
       do k=1,n2-1
        ar=s(k)*cos(2.*pi*r(k))
        ai=s(k)*sin(2.*pi*r(k))
        lam(k)=cmplx(ar,ai)
        lam(mp-k)=cmplx(ar,-ai)
       enddo
        else
         do k=1,n2
          ar=s(k)*cos(2.*pi*r(k))
          ai=s(k)*sin(2.*pi*r(k))
          lam(k)=cmplx(ar,ai)
          lam(mp-k)=cmplx(ar,-ai)
         enddo
      endif
c----------
c  System matrix
      do k=0,mp-1
       z(k)=lam(k)
      enddo
      call sfftf(mp,z,v)
      do k1=1,m
       do k2=1,mp
        if(k1>=k2) then
         b(k1,k2)=real(z(k1-k2))/float(mp)
          else
           b(k1,k2)=real(z(mp+k1-k2))/float(mp)
        endif
       enddo
      enddo
c-------
      init='n'
c  Start recursion
      do t=1,p
       call l90rand(m,e,ch,seed,init,io)
       do km=1,m
        i=(km-1)*n+t
        x(i)=sig(km)*e(km)
       enddo
      enddo
c  Rest of sample
      do t=p+1,n
       call l90rand(m,e,ch,seed,init,io)
c  Shift x(t-1),...,x(t-p) to s
       do j=1,p
        do km=1,m
         i=(km-1)*n+t-j
         s((j-1)*m+km)=x(i)
        enddo
       enddo
       call markov(s,lam,z,v,t,mp,io)
       do km=1,m
        i=(km-1)*n+t
        x(i)=s(km)+sig(km)*e(km)
       enddo
      enddo
c----------
      call ahead(x,y,b,e,sig,s,ch,m,p,q,n,seed,io)
c  Call varm
      if(m<=40) then
       call varm(x,c,s,r,r2,m,p,n,ier,io)
       if(ier>0) then
        write(io,*) 'Error in subroutine varm ier = ',ier
        stop
       endif
        else
         call varsimple(x,c,s,m,p,n,ier,io)
         if(ier>0) then
          write(io,*) 'Error in subroutine varsimple ier = ',ier
          stop
         endif
      endif
c----------
c  Write out VAR coefficients
      do k1=1,m
       do kp=1,p
        do k2=1,m
         jp=(kp-1)*m+k2
         tst=sngl(c(jp,k1))/s(jp)
         if(tst>1.e4) then
          tst=1.e4
           elseif(tst<(-1.e4)) then
            tst=-1.e4
         endif
         write(io,'(2x,''c('',i2,'','',i2,'') = '',f8.3,2x,''c^ = '',f8.
     *3,2x,''se(a) = '',f8.4,2x,''t = '',f8.2)') k1,jp,b(k1,jp),c(k1,jp)
     *,s(jp),tst
        enddo
       enddo
      enddo
c----------
      do k=1,m
       call stat(n,x((k-1)*n+1),am,sd,sk,c4,c6,xmx,xmn,io)
       e(k)=sd
       if(m<=40) then
        write(io,'(50(''=''))')
        write(io,'(7x,''Statistics for x('',i2,'')''/)') k
        write(io, '(7x,''Mean = '',g9.2,2x,''Sig = '',g12.4)') am,sd
        write(io,'(7x,''Skewness = '',g12.3,2x,''Kurtosis = '',f12.2)')
     *sk,c4
        write(io,'(7x,''Max ='',g17.5,2x,''Min ='',g17.5)') xmx,xmn
        write(io,'(50(''*''))')
c-----------
        call stat(n,r((k-1)*n+1),rm,rs,rw,r4,r6,rmx,rmn,io)
        write(io,'(7x,''Statistics for x('',i2,'') residuals''/)') k
c  Write R square
        write(io,'(7x,''Adjusted R Squared = '',f7.3/)') r2(k)
        write(io, '(7x,''Mean = '',g9.2,2x,''Sig = '',f9.4,2x,''True Sig
     * = '',f9.4)') rm,rs,sig(k)
        write(io,'(7x,''Skewness = '',g12.3,2x,''Kurtosis = '',f12.2)')
     *rw,r4
        write(io,'(7x,''Max ='',g17.5,2x,''Min ='',g17.5
     *)') rmx,rmn
       endif
      enddo
c---------
      if(tr=='t') then
c  Use true VAR coefficients for forecasting
       do k1=1,m
        do kp=1,p
         do k2=1,m
          jp=(kp-1)*m+k2
          c(k1,jp)=dble(b(k1,jp))
         enddo
        enddo
       enddo
      endif
      call forvar(x,r,c,s,m,p,q,n,io)
c  Write out forecasts errors
      write(io,'(50(''*'')/)')
      write(io,'(7x,''Forecast Errors Relative to the Standard Deviation
     * of the Equation '',i3,'' steps ahead''/)') q
      write(io,'(''Eqn # '')')
      write(io,*)
      do k=1,m
       write(io,'(i3,:,99(1x,f7.2))') k,((r((q-t)*m+k)-y((q-t)*m+k))/e(k
     *)*100.,t=1,q)
      enddo
      write(io,'(3x,:,99(3x,i3,2x))') (t,t=1,q)
      write(io,'(37(''=''))')
c-------
      if(imp=='i') then
       call pulse(y,b,r,m,p,q,n,7)
       do t=1,q
        write(7,'(1x,:,99(f7.3,1x))') (y((q-t)*m+k),k=1,m)
       enddo
       write(7,'(/:,99(3x,i3,2x))') (k,k=1,m)
       write(7,'(/7x,''Impulse Responses of Variates'')')
      endif
      return
c----------
      end
c
c===========================================
c
      subroutine varm(x,c,s,r,r2,m,p,n,ier,io)
c********************************************
c  Compute coefficients of var with m variables and p lags
c  M. J. Hinich  Version 11-8-99
c  Index i for the ith equation with m*p variates
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+ cim2*xm(t-2)
c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Iput: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c  Output: c - VAR coefficients in row major m*m*p array
c   s - standard errors of c coeff's
c   r - residuals for the m variables in column major order
c   r2 - m adjusted R squares
c   ier = 1: A is not positive definite, ier = 2: m < 1,
c   ier = 3: p < 1, ier = 4: m*m*p > n/2
c============================================
      allocatable::a,sg
      integer p,mar
      real x(m*n),r(m*n),s(m*m*p),r2(m),sg(:)
      real*8 c(m*m*p),a(:)
      intent(in) x,m,p,n,io
      intent(out) c,r,r2,ier
c********************************************
c  Traps
      if(m<1) then
       write(io,*)'m ',m,' < 1'
       ier=2
       return
        elseif(p<1) then
         write(io,*)'p ',p,' < 1'
         ier=3
         return
          elseif(m*m*p>n/2) then
           write(io,*) 'm*m*p > n/2 = ',n/2,' m = ',m,' p = ',p
           ier=4
           return
      endif
c----------
      nc=m*m*p
      i2=nc*nc
      i3=i2+nc
      ia=i3+nc
c----------
c  Allocate space
      allocate(a(0:ia-1),sg(0:m-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in varm'
       stop
      endif
c  Initalize
      do k=1,ia-1
       a(k)=0.d0
      enddo
      call var(x,r,a(0),a(i2),a(i3),sg,m,p,n,nc,c,s,r2,ier,io)
c----------
      deallocate(a,sg)
      return
      end
c
c===========================================
c
      subroutine var(x,r,a,b,d,sg,m,p,n,nc,c,s,r2,ier,io)
c********************************************
c  Least squares fit of a VARm(p) model. Version 8-19-97
c  Iput: x - data matrix of n observations
c  Output: r - residuals, c - VAR coefficients
c   s - standard errors of coefficients
c============================================
      integer p,t
      real x(m*n),r(m*n),s(nc),r2(m),sg(m)
      real*8 a(nc,nc),b(nc),d(nc),c(nc),sm,dn,a2,st
      intent(in) x,m,p,n,nc,io
      intent(out) c,r,r2,ier
c********************************************
      dn=dfloat(n)
      mp=m*p
c  Covariance of xk1(t) & xk2(t-kp)
      do k1=1,m
       do kp=0,p
        do k2=1,m
c  Average over time
         sm=0.d0
         do t=1,n-kp
          sm=sm+dble(x((k1-1)*n+t+kp))*dble(x((k2-1)*n+t))
          if(kp==0 .and. k1==k2) then
c  Sum of squares of xk1(t)
           r2(k1)=sngl(sm)
          endif
         enddo
         sm=sm/dn
c  b array
         if(kp>0) then
          jp=(k1-1)*mp+(kp-1)*m+k2
          b(jp)=sm
         endif
c  Reproduce sm to fill up the matrix a
         do jm=1,m
          if(kp<p) then
           do jp=0,p-1-kp
            ip=(jm-1)*mp
            j1=jp*m+ip+k1
            j2=(jp+kp)*m+ip+k2
            a(j1,j2)=sm
           enddo
          endif
         enddo
        enddo
       enddo
      enddo
c----------
c  Cholesky decomposition of a
      call dcholdc(a,nc,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       ier=1
       return
      endif
c----------
      call dcholsl(a,nc,d,b,c)
c  Compute residuals & standard deviations
      do k1=1,m
      a2=0.d0
c  Zero out first p residuals
       do t=1,p
        r((k1-1)*n+t)=0.
       enddo
c----------
       do t=p+1,n
        sm=0.d0
c  Residuals
        do kp=1,p
         do k2=1,m
          jp=(k1-1)*mp+(kp-1)*m+k2
          kx=(k2-1)*n+t-kp
          sm=sm+c(jp)*dble(x(kx))
         enddo
        enddo
        st=dble(x((k1-1)*n+t))-sm
        r((k1-1)*n+t)=sngl(st)
        a2=a2+st*st
       enddo
c  End of loop on t
       nap=n-mp-1
       sg(k1)=sngl(dsqrt(a2/dfloat(nap)))
c  Adjusted R squares
       rsq=(sngl(a2)/r2(k1))/(1.-float(mp)/float(n-1))
       r2(k1)=1-rsq
      enddo
c  End of loop on k1
c  Inverse of A
      call dinverse(a,nc,d)
c  Standard errors of c's
      do k1=1,m
       do kp=1,p
        do k2=1,m
         jp=(k1-1)*mp+(kp-1)*m+k2
         s(jp)=sg(k1)*sngl(dsqrt(a(jp,jp)))
        enddo
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine varsimple(x,c,s,m,p,n,ier,io)
c********************************************
c  Compute coefficients of var with m variables and p lags
c  M. J. Hinich  Version 1-19-00
c  Index i for the ith equation with m*p variates
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+ cim2*xm(t-2)
c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Iput: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c  Output: c - VAR coefficients in row major m*m*p array
c   s - standard errors of c coeff's
c   ier = 1: A is not positive definite, ier = 2: m < 1,
c   ier = 3: p < 1, ier = 4: m*m*p > n/2
c============================================
      allocatable::a,sg
      integer p,mar
      real x(m*n),s(m*m*p),sg(:)
      real*8 c(m*m*p),a(:)
      intent(in) x,m,p,n,io
      intent(out) c,ier
c********************************************
c  Traps
      if(m<1) then
       write(io,*) 'm ',m,' < 1'
       ier=2
       return
        elseif(p<1) then
         write(io,*) 'p ',p,' < 1'
         ier=3
         return
          elseif(m*m*p>n/2) then
           write(io,*) 'm*m*p > n/2 = ',n/2,' m = ',m,' p = ',p
           ier=4
           return
      endif
c----------
      nc=m*m*p
      i2=nc*nc
      i3=i2+nc
      ia=i3+nc
c----------
c  Allocate space
      allocate(a(0:ia-1),sg(0:m-1),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space in varsimple'
       stop
      endif
c  Initalize
      do k=1,ia-1
       a(k)=0.d0
      enddo
      call varx(x,a(0),a(i2),a(i3),sg,m,p,n,nc,c,s,ier,io)
c----------
      deallocate(a,sg)
      return
      end
c
c===========================================
c
      subroutine varx(x,a,b,d,sg,m,p,n,nc,c,s,ier,io)
c********************************************
c  Least squares fit of a VARm(p) model. Version 1-19-00
c  Iput: x - data matrix of n observations
c  Output: c - VAR coefficients, s - standard errors of coefficients
c============================================
      integer p,t
      real x(m*n),s(nc),sg(m)
      real*8 a(nc,nc),b(nc),d(nc),c(nc),sm,dn,a2,st
      intent(in) x,m,p,n,nc,io
      intent(out) c,ier
c********************************************
      dn=dfloat(n)
      mp=m*p
c  Covariance of xk1(t) & xk2(t-kp)
      do k1=1,m
       do kp=0,p
        do k2=1,m
c  Average over time
         sm=0.d0
         do t=1,n-kp
          sm=sm+dble(x((k1-1)*n+t+kp))*dble(x((k2-1)*n+t))
         enddo
         sm=sm/dn
c  b array
         if(kp>0) then
          jp=(k1-1)*mp+(kp-1)*m+k2
          b(jp)=sm
         endif
c  Reproduce sm to fill up the matrix a
         do jm=1,m
          if(kp<p) then
           do jp=0,p-1-kp
            ip=(jm-1)*mp
            j1=jp*m+ip+k1
            j2=(jp+kp)*m+ip+k2
            a(j1,j2)=sm
           enddo
          endif
         enddo
        enddo
       enddo
      enddo
c----------
c  Cholesky decomposition of a
      call dcholdc(a,nc,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       ier=1
       return
      endif
c----------
      call dcholsl(a,nc,d,b,c)
c  Compute residuals & standard deviations
      do k1=1,m
      a2=0.d0
       do t=p+1,n
        sm=0.d0
c  Residuals
        do kp=1,p
         do k2=1,m
          jp=(k1-1)*mp+(kp-1)*m+k2
          kx=(k2-1)*n+t-kp
          sm=sm+c(jp)*dble(x(kx))
         enddo
        enddo
        st=dble(x((k1-1)*n+t))-sm
        a2=a2+st*st
       enddo
c  End of loop on t
       nap=n-mp-1
       sg(k1)=sngl(dsqrt(a2/dfloat(nap)))
      enddo
c  End of loop on k1
c  Inverse of A
      call dinverse(a,nc,d)
c  Standard errors of c's
      do k1=1,m
       do kp=1,p
        do k2=1,m
         jp=(k1-1)*mp+(kp-1)*m+k2
         s(jp)=sg(k1)*sngl(dsqrt(a(jp,jp)))
        enddo
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine forvar(x,y,c,w,m,p,q,n,io)
c********************************************
c  Forecasts q steps ahead for a VAR(m,p) model
c  M. J. Hinich  Version 8-9-97
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+
c   cim2*xm(t-2) +...+ c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Iput: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c   c - VAR coefficients in m*m*p array.
c  Output: y - x1(i),...,xm(i) - i point forecasts, i = q,q-1,..,1
c============================================
      integer m,p,q,n,t,tc
      real x(m*n),y(m*q),w(m)
      real*8 c(m,m*p),sm
      intent(in) x,c,m,p,q,n,io
      intent(out) y
c********************************************
      do t=1,m*q
       y(t)=0.0
      enddo
c Shift the last p observations to y
      do t=n,n-p+1,-1
       do k=1,m
        i=(n-t)*m+k
        y(i)=x((k-1)*n+t)
       enddo
      enddo
c----------
      do t=1,q
       do k1=1,m
        sm=0.d0
        do kp=1,p
         do k2=1,m
          kx=(kp-1)*m+k2
          sm=sm+c(k1,kx)*dble(y(kx))
         enddo
        enddo
c kc forecast
        w(k1)=sngl(sm)
       enddo
c Foreshift the y's
       do ts=q-1,1,-1
        do k1=1,m
         kt=(ts-1)*m+k1
         y(kt+m)=y(kt)
        enddo
       enddo
c Shift kc forecast to first block
       do k=1,m
        y(k)=w(k)
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine ahead(x,y,b,e,sig,w,ch,m,p,q,n,seed,io)
c********************************************
c  VAR q steps ahead
c  M. J. Hinich  Version 11-9-99
c  xi(t)=ci11*x1(t-1) +...+ cim1*xm(t-1) + ci12*x1(t-2) +...+ cim2*xm(t-
c  + ci1p*x1(t-p) +...+ cimp*xm(t-p) + ei(t)
c  Iput: x - data: x1(1),...,x1(n),x2(1),...,x2(n),...,xm(1),...,xm(n)
c   b - VAR coefficients
c  Output: y - x1(i),...,xm(i) - i steps ahead, i = 1,...,q
c============================================
      integer m,p,q,n,seed,t,tc
      real x(m*n),y(m*q),b(m,m*p),e(m*q),sig(m),w(m)
      character ch,init
      intent(in) x,b,sig,ch,m,p,q,n,seed,io
      intent(out) y
c********************************************
      init='n'
      call l90rand(m*q,e,ch,seed,init,io)
c  Shift the last p observations to y
      do t=n,n-p+1,-1
       do k=1,m
        i=(n-t)*m+k
        y(i)=x((k-1)*n+t)
       enddo
      enddo
c----------
      do t=1,q
       do k1=1,m
        ip=(t-1)*m+k1
        sm=0.
        do kp=1,p
         do k2=1,m
          jp=(kp-1)*m+k2
          sm=sm+b(k1,jp)*y(jp)
         enddo
        enddo
        w(k1)=sm+sig(k1)*e(ip)
       enddo
c  Foreshift the y's
       do tc=q-1,1,-1
        do k1=1,m
         kt=(tc-1)*m+k1
         y(kt+m)=y(kt)
        enddo
       enddo
       do k=1,m
        y(k)=w(k)
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine pulse(y,b,w,m,p,q,n,io)
c********************************************
c  Impulse response of VAR for q steps ahead
c  M. J. Hinich  Version 5-27-97
c  Output: y - x1(i),...,xm(i) - i steps ahead, i = 1,...,q
c============================================
      integer m,p,q,n,t,tc
      real y(m*q),b(m,m*p),w(m)
      intent(in) b,m,p,q,n,io
      intent(out) y
c********************************************
      do k=1,m*q
       y(k)=0.
      enddo
c  Unit impusle for each x
      do k=1,m
       y(k)=1.
      enddo
      do t=1,q
       do k1=1,m
        sm=0.
        do kp=1,p
         do k2=1,m
          jp=(kp-1)*m+k2
          sm=sm+b(k1,jp)*y(jp)
         enddo
        enddo
        w(k1)=sm
       enddo
c  Foreshift the y's
       do tc=q-1,1,-1
        do k1=1,m
         kt=(tc-1)*m+k1
         y(kt+m)=y(kt)
        enddo
       enddo
       do k=1,m
        y(k)=w(k)
       enddo
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine markov(x,lam,z,v,t,mp,io)
c********************************************
c  Linear system using the rotation matrix t=(sin(2*pi*jk/n) on the
c  uncoupled complex AR(1) solutions
c  M. J. Hinich  Version 1-14-00
c============================================
      integer t,mp,io
      real x(mp)
      complex z(0:mp-1),lam(0:mp-1),v(4*mp+15)
      intent(in) lam,v,t,mp,io
      intent(out) x
c********************************************
      do k=1,mp
       z(k-1)=cmplx(x(k)/float(mp),0.)
      enddo
      call sfftf(mp,z,v)
c  Contraction by lam & conjugate
      do k=0,mp-1
       z(k)=conjg(lam(k)*z(k))
      enddo
      call sfftf(mp,z,v)
      do k=1,mp
       x(k)=real(z(k-1))
      enddo
c----------
      return
      end
c
c===========================================
c
      subroutine stat(n,x,mean,sd,sk,c4,c6,max,min,io)
c**************************************************
c     M.J. Hinich and J.W. Dalle Molle
c     Version : 8-15-97
c     Precision : single precision/double accumulation
c
c     Function : Calculates mean, standard deviation, skewness - sk,
c      kurtosis - c4, 6-th order cumulant - c6, maximum - max, and
c      minimum - min of the sample.
c
c  Input: n - sample size, x - data array
c
c=================================================
      real*8 s1,s2,s3,s4,s6,dx,dx3,dn,xmin,xmax
      real x(n),mean,max,min
      intent(in) n,x,io
      intent(out) mean,sd,sk,c4,c6,max,min
c**************************************************
c Initialize accumulators
      xmax=x(1)
      xmin=x(1)
      s1=0.d0
      s2=0.d0
      s3=0.d0
      s4=0.d0
      s6=0.d0
      dn=dfloat(n)
      do k=1,n
       dx=dble(x(k))
       xmax=dmax1(dx,xmax)
       xmin=dmin1(dx,xmin)
       s1=s1+dx
       s2=s2+dx*dx
      enddo
c----------
      max=xmax
      min=xmin
      s1=s1/dn
      mean=s1
c------------
      sd=sngl(s2)-float(n)*mean*mean
c  Check if variance = 0
      if(sd<=0.) then
       write(io,*) 'Variance in STAT = 0 & is set = 1.'
       sd=1.
       return
      endif
c  Compute standard deviation
      rn=float(n-1)
      sd=sqrt(sd/rn)
c---------------------
c  Compute other stats
      do k=1,n
       dx=dble(x(k))-s1
       dx3=dx*dx*dx
       s3=s3+dx3
       s4=s4+dx3*dx
       s6=s6+dx3*dx3
      enddo
c----------
      dx3=sd*sd*sd
      sk=s3/(dn*dx3)
      c4=s4/(dn*dx3*sd)-3.d0
      s6=s6/(dx3*dx3)
      c6=s6/dn-15.d0*c4-10.d0*sk*sk-15.d0
c-------
      return
      end