      PROGRAM MTARFOR
      USE numerical_libraries
cc      
cc  FORTRAN77 code
cc  Coded by Ruey S. Tsay and modified by JDG (2007)
cc
cc  one-step ahead forecasting from bivariate TAR model
cc  program consists of 3 parts:
cc
cc  1) mtarch.for (subroutine):   finding optimal threshold value
cc  2) mtares.for (main program): estimating parameters
cc  3) one-step ahead forecasting + computing smooth tests
cc  
cc  Reference:
cc  De Gooijer, J.G. (2007).
cc     Power of the Neyman smooth test for evaluating multivariate 
cc     forecast densities.
cc     Journal of Applied Statistics, 34(4), 371-382.
cc     DOI: 10.1080/02664760701231526.
cc***********************************************************
cc  This program is intended for nonlinear time series analysis.
cc  It estimates a multivariate TAR model.
cc
cc  Created on May 20, 1997 to analyze intraday data sets.
cc  Modified on July 7 so that each regime may have its own AR order.
cc
cc  Further modified on November 14, 1997 to allow for diferent component 
cc  model in each regime. The model is now very flexible.
cc
cc  The program estimates a threshold-regression model: threshold variable
cc  and its threshold must be given.
cc  input: dat(.,.) is the data matrix
cc    ,,   dt(.,.)  is original data matrix
cc  locy: the column number in dat(.,.) that gives the dependent var.
cc  locthr: column in which threshold variable is stored in dat(.,.)
cc  locx: columns where independent variables are located in dat(.,.).
cc  iorx: order of each independent variables
cc  ixlag: lags of each independent variables
cc  p: order of the dependent variable (fixed at 8 for all regimes)
cc  nx: number of exogeneous variables.
cc  inx: number of regressors due to exogeneous variables.
cc  (nst, nend): data span used in the program
cc         
cccccccccccccc Current settings   *********************
c    max-variables = 5
cc   max-ar-p = 29
cc   max-obs = 10000
cc   max-number-parameters for each component = 50
cc
cc   Maximum number of regimes = 5
c
       parameter(maxn=10000,md=50,mk=5)
       real*8 xpx(md,md),xpy(md),xpxinv(md,md) 
       real*8 phi(md),y(maxn),dat(maxn,6), thr(5)
       real*8 resi(maxn,mk), res(maxn,mk), tem, chk
       real*8 x(maxn,md), rms(mk), ratio,std,gthr(5)
       real*8 aic, bic, dlog, wk(md,md), rrcov(mk,mk)
       real*8 stdres(maxn,mk), sigma(mk), tmp,ipos2
       real*8 phia(3,2,20),pred(5),resia(3,maxn,2),rnn(2)
       real*8 predtot(maxn,mk),z(mk),sigm(mk,mk)
       real pcdf(maxn),cdfr(maxn)
       real xx1,xx2,zmar,zjoint,uu1(5),uu2(5),uu3,uu4,sprob(5)
       real rho1,rho2,rho11,ymar(maxn),yjoint(maxn)
       real rhon1(mk),rhon11(mk),rhon2(mk)
       real mu1,mu2,mu3,mu4,pr(1000,2)
       real cdfa,cdfb,cdff
       real dens1(1000),dens2(1000)
       real dt(maxn,6)
c
       integer nob, p, ip, lagth, nreg, npts(5), inx, n3
       integer locy(5), locthr, nx, ny, locx(5), nvar, i, j 
       integer iorx(5), ist, ia, ii, jj, it
       integer icnt, nst, nend, idx, iar(5)
       integer ixlag(15,5), morx, icom,seed,ipos,s
c
       character nfile*30
c
       iout = 17
       seed=1234567
c
c---------- nob = # of observations; p = ar-order
c
c       write(6,1)
c    1  format(1x,'nob, nvar, n_y, n_x, delay, nregime, loc_thr: ',$)
c       read(5,*)nob, nvar, ny, nx, lagth, nreg, locthr

c***************   fixed values   ****************
       n3=7060     ! total length of series
       nvar=3
       ny=2
       nx=1
       lagth=1  ! delay
       nreg=2   ! number of regimes
       locthr=3 ! LOC THRESHOLD VAR in data file (variable X, col 3) 
c ***********************************************
	 
       if(n3 .le. 0)stop
       if(n3 .gt. maxn)n3 = maxn
       if(nvar .le. 0)nvar = 1
       if((ny+nx).gt.nvar)then
        print*,'Number of variables fails to match'
        stop
       endif
c       if(nvax .gt. 5)then
c        print*,'too many variables! (the maximum is 5)'
c        stop
c       endif
       if(lagth.lt.0)lagth = 0
       if(locthr.le.0)locthr = 1
c
       write(6,2)
    2  format(1x,'Data_file Name: ',$)
       read(5,3)nfile
    3  format(a30)
       open(unit=21,file=nfile,status='old')
       do 4 i = 1, n3      
 4        read(21,*)(dt(i,j),j=1,nvar)
       close(21)
       
c       do 1 i=1,n3
c1      write(17,*) (dt(i,j),j=1,nvar)
c
c       write(6,401)
c 401   format(1x,'input data span (1 = yes): ',$)
c       read(5,*) ii
        ii=0
        if(ii.eq.1)then
        write(6,402)
        read(5,*) nst, nend
        if(nend.gt.n3) nend=n3
        if(nst.lt.1) nst = 1
        n3 = nend-nst+1
        if(nst.eq.1)go to 415
        do 410 it=1, n3
         do 408 i=1, nvar
         write(17,*) ' dt ', dt(it-1+nst,i),it,i,nst
 408        dat(it,i) = dt(it-1+nst,i)
 410        continue
       endif
 402   format(1x,'data span (nst, nend): ',$)
 415   continue
c
c      write(6,5)
c 5    format(1x,'loc_y: ',$)
c      read(5,*) (locy(i),i=1,ny)
       locy(1)=1
       locy(2)=2
       locx(1)=3
c
c      if(nx.gt.0)then
c       write(6,7)
c       read(5,*)(locx(i),i=1,nx)
c      endif
c 7      format(1x,'loc_x: ',$)
c
c *********  change the next line  ***********************
c
c      s=100  ! 100  !  or 300
	s=200
c ********************************************************

      WRITE(17,*) 'PROGRAM: MTARFOR.FOR in map TSAY FORTRAN'
      WRITE(17,*) 'NOB = ', N3     ! LENGTH OF SERIES
      WRITE(17,*) 'NVAR = ', NVAR ,'NY = ', NY, 'NX = ',nx 
      WRITE(17,*) 'DELAY (LAGTH) = ', LAGTH 
      WRITE(17,*) 'NREGIME = ', NREG
      WRITE(17,*) 'LOCTHR =', LOCTHR 
      WRITE(17,*) 'S = ', S
      WRITE(17,*) 

      iprob1=0
      iprob2=0
      
      do 99999 iwq=1,s    ! big loop for rolling forecasts
      
      nob=n3-(s+1)+iwq ! N3 = 7060, (s+1)=LAST IN-SAMPLE
        
   !  FOR IWQ=1, Nob=7060-101+1=6960 ESTIMATION:
c     with s=2 iwq=7060-2=7058 last in sample value; forec 7059
c     with s=2 iwq=7060-1=7059 last in sample; forecast 7060

   !  FOR IWQ=2, Nob=7060-101+2=6961 ESTIMATION
   !  ...
   !  FOR IWQ=100, Nob=7060-101+100=7059 
   
       write(*,*) 'iwq (loop) =', iwq                         
       do 193 it=1,nob             
            do 199 i=1,nvar
c            write(*,*) dt(it,i),it,i
199            dat(it,i)=dt(it,i) 
193         continue  
       write(*,*) ' dat ', dt(nob+1,1),dt(nob+1,2),nob+1

       if(nreg.gt.1)then 
c--------- input thresholds
cc         write(6,11) nreg-1
cc         read(5,*)(thr(i),i = 2, nreg)
       call mtarsch(nob,dat,gthr)
       thr(2)=gthr(1)
c       write(*,*) 'thr ', thr(2)
       
       endif
        thr(1) = -9.9d10
        i = nreg+1
        thr(i) = 9.9d10
   11  format(1x,'input', i5,' thresholds: ',$)
c
c       write(6,8)
c 8     format(1x,'max-ar-order & max-exog-order: ',$)
c       read(5,*) p, morx
       p=8
       morx=1
       ist = max0(p,morx,lagth)+1
c
c------In what follows, perform estimation eqn by eqn & 
c------regime by regime.
c
       aic = 0.0d0
       bic = 0.0d0
c
      do 8000 ia=1, nreg
c
       do 7000 icom = 1, ny
c
cc       print*,'regime: ',ia,' component of focus: ',icom
cc       write(6,61) 
cc        read(5,*) (iar(it),it=1,ny)
c61     format(1x,'AR order by components: ',$)
       iar(1)=8
       iar(2)=8
c
       if(nx.gt.0)then
c       write(6,9)
c       read(5,*)(iorx(i),i=1,nx)
       iorx(1)=1
c------------------------- give the default lags of inde. vars.
        do 44 i=1, nx
         if(iorx(i).gt.15)then
          print*,'max[#(lags)] is 15. Reset to 15.'
          iorx(i) = 15
         endif
c
         do 42 j=1, iorx(i)
 42         ixlag(j,i) = j
 44         continue
c
c       write(6,56)
c       read(5,*) it
        it=0
        
        if(it.eq.1)then
         do 54 i=1, nx
          write(6,57) i
          read(5,*)(ixlag(j,i),j=1,iorx(i))
          do 53 j=1, iorx(i)
           if(ixlag(j,i).lt.0)then
              print*,'negative order not allowed'
              stop
             endif
 53        continue
 54       continue
        endif
       endif
 9      format(1x,'#(lags) of the x-variables: ',$)
 56     format(1x,'input detailed lags of exog. vars? (1=y): ',$)
 57     format(1x,'lags of',i3,'-th exo. var: ',$)
c      
c----------------- inx = # of regressors from x's variables.
      inx = 0
      if(nx.gt.0)then
       do 17 i=1, nx
          inx = inx + iorx(i)
 17       continue
      endif
c
c------- including a constant term (ip = dimension of x-matrix)
       ip = 1+ inx
       do 18 i=1, ny
 18       ip = ip + iar(i)
       if(ip.gt.md)then
        print*,'dimension exceeds program specification'
        stop
       endif
c------------------- icnt: number of data points in this regime.
      icnt = 0
        do 100 it = ist, nob
         chk = dat(it-lagth,locthr) 
c         write(*,*) 'thr', thr(ia), ia
         
         if(chk.lt.thr(ia))go to 100
         if(chk.ge.thr(ia+1))go to 100
c------------- the observation is in the ia-th regime
         icnt = icnt+1
c         write(*,*) ' ist ', ist,nob,chk,it,icnt         
         y(icnt) = dat(it,locy(icom))
         x(icnt,1) = 1.0d0
         idx = 1
         do 32 jj=1, ny
          if(iar(jj).le.0)go to 32
           do 30 i=1, iar(jj)
             idx=idx+1
 30          x(icnt,idx) = dat(it-i,locy(jj))
 32        continue
c
        if(nx.gt.0)then
         do 50 ii=1, nx
          if(iorx(ii).gt.0)then
           do 40 jj=1, iorx(ii)
              idx=idx+1
 40           x(icnt,idx)=dat(it-ixlag(jj,ii),locx(ii))
          endif
 50      continue
        endif
c
 100    continue
c
        npts(ia) = icnt
c        write(*,*) 'icnt = ', icnt
        if(icnt.le.ip)then
         print*,'insufficient data points in the', ia,'-th regime'
         stop
        endif
c
       if(idx.ne.ip)then
        print*,'dimension error'
        stop
       endif
       call regr(y,x,maxn,md,ip,icnt,xpx,xpxinv,xpy,phi)
c
c----------residuals & residual variance for component ''icom''
       tmp = 0.0d0
       do 120 it=1, icnt  ! 1,...,# data points in regime
         tem = 0.0d0
         do 110 ii=1, ip
 110      tem = tem + phi(ii)*x(it,ii)
 
       tem = y(it)-tem
       resi(it,icom) = tem
       resia(ia,it,icom)=tem
c       write(iout,*) resia(ia,it,icom),ia,it,icom
       tmp = tmp + tem*tem
 120   continue     
c
c       write(iout,152) ia,npts(ia)
c
       rms(icom) = tmp/dfloat(icnt-ip)
       sigma(icom) = dsqrt(rms(icom))
c
c        write(iout,301) icom
        do 150 ii=1, ip
         std = dsqrt(xpxinv(ii,ii)*rms(icom))
         ratio = phi(ii)/std
c        write(iout,151) phi(ii), std, ratio 
 150    continue
 300    continue
c 151    format(1x,4f10.5)
c 301    format(1x,'component ', i3)
c
c 152    format(1x,'Regime:',i4,' Nobs:',i5)
c
          do 1553 ii=1,ip
          phia(ia,icom,ii)=phi(ii)
c          write(iout,*) phia(ia,icom,ii),ia,icom,ii
1553       continue
        
 7000   continue  ! end loop icom=1,ny
c----------------- Residual covariance matrix (MLE-type and LS-type)
       do 200 i=1, ny
        do 190 j=1, i
         tem = 0.0d0
         do 180 it=1, icnt
c         write(*,*) resi(it,j),it,j
 180        tem=tem+resi(it,i)*resi(it,j)
        rrcov(i,j) = tem/dfloat(icnt-ip)
        rrcov(j,i) = rrcov(i,j)
c        write(*,*) 'rrcov', rrcov(i,j),i,j,icnt
        xpx(j,i) = tem/dfloat(icnt)
 190    xpx(i,j) = tem
 200    continue
c
       call detmtx(xpx,wk,md,ny,tem)
**
       tem = dlog(tem)*dfloat(icnt)
c       write(iout,201) tem
       aic = aic + tem + 2.0d0*dfloat(ip*ny)
       bic = bic + tem + dlog(dfloat(icnt))*dfloat(ip*ny)
c
c 201   format(1x,'(log(det)*nobs) = ',f16.6)
c
c        write(iout,302)
c        do 160 i=1, ny
c  160      write(iout,153)(rrcov(i,j),j=1,ny)
c  153      format(1x,5f20.10)
c 302     format(1x,'Residual cov-mtx: ')
c
      sigm(ia,1)=sqrt(rrcov(1,1))
      sigm(ia,2)=sqrt(rrcov(2,2))
      rhon1(ia)=rrcov(1,2)/(sigm(ia,1)*sigm(ia,2)) ! rho
      rhon2(ia)=sqrt(1-rhon1(ia)*rhon1(ia))
      rhon11(ia)=sigm(ia,1)*rhon2(ia)
c     write(17,*) sigm(ia,1),sigm(ia,2),rhon1(ia),ia
c
c------ store residuals and standardized residuals
        icnt = 0
        do 500 it=ist, nob
         chk = dat(it-lagth,locthr)  
         if(chk.lt.thr(ia))go to 500
         if(chk.ge.thr(ia+1))go to 500
         icnt = icnt+1
         do 450 j=1, ny
            stdres(it,j) = resi(icnt,j)/sigma(j)
 450        res(it,j) = resi(icnt,j)
 500     continue
c
c--------------- return for another regime
 8000  continue  ! end loop ia=1,nreg
c
c       write(iout,3002)(thr(i),i=2,nreg)
c 3002  format(1x,'Thresholds: ',5f12.6)
c       write(iout,3003)aic,bic
c       write(6,3003)aic,bic
c 3003  format(1x,'Overall AIC & BIC',2e16.8)
c

c
c       do 3100 it=ist, nob
c          write(18,3101) (stdres(it,j),j=1,ny)
c 3100     write(16,3101) (res(it,j),j=1,ny)
c 3101     format(1x,5f10.6)
c
c       print*,'residuals in fort.16; std. residuals in fort.18'

c      ******* one-step ahead forecasting **********
c      so it=nob+1 with lagth=1 gives (it-lagth)=nob

       chk=dat(nob,locthr)
c       write(iout,*) ' chk (threshold value)= ',chk, thr(2)
       
       if(chk.lt.thr(2)) then  ! switch betwn regimes: < thr(2)
       tem=0.0d0
       do 444 icom=1,ny  ! ny=2 
       
         do 443 ii=1,8   ! AR order each regime fixed at 8(max)
         tem=tem+phia(1,icom,2*ii)*dat(nob-ii+1,1)+
     *  phia(1,icom,2*ii+1)*dat(nob-ii+1,2)
c         write(iout,*) dat(nob-ii+1,1),dat(nob-ii+1,2),dat(nob,3),nob,ii
443      continue  ! end ii=1,8

         pred(icom)=tem+phia(1,icom,1)+phia(1,icom,18)*dat(nob,3)
c        + constant term
c       write(iout,*) ' pred1 = ', pred(icom),icom      
444    continue           ! end loop icom=1,ny
       sigma(1)=sigm(1,1)
       sigma(2)=sigm(1,2)
       rho11=rhon11(1)
       rho1=rhon1(1)      ! rho
       
       else    ! switch between regimes, in this case > thr(2)
       
       tem=0.0d0
       do 446 icom=1,ny  ! ny=2 two variables
         do 447 ii=1,8   ! AR-order is fixed at 8
         tem=tem+phia(2,icom,2*ii)*dat(nob-ii+1,1)+
     *   phia(2,icom,2*ii+1)*dat(nob-ii+1,2)
      
c        write(iout,*) dat(nob-ii+1,1),dat(nob-ii+1,2),nob,ii,dat(nob,1),dat(nob,2)
447      continue  ! end loop ii=1,8
      pred(icom)=tem+phia(2,icom,1)+phia(2,icom,18)*dat(nob,3)
c      write(iout,*) ' pred2 = ', pred(icom),icom
446   continue     ! end loop icom=1,ny
      sigma(1)=sigm(2,1)
      sigma(2)=sigm(2,2)
      rho11=rhon11(2) ! (1-rho^2)*sigma^2
      rho1=rhon1(2)   ! rho 
      endif           ! end switch between thresholds
      
c      write(17,4444) pred(1),pred(2),dt(nob+1,1),dt(nob+1,2)
c4444  format(4f20.10)

      pr(iwq,1)=pred(1)
      pr(iwq,2)=pred(2)
c
c     **** standardize predictions ********

c      write(*,*) ' rho ', rho1
c      write(17,*) 'rho , sigma1 and sigma2 ='
c      write(17,881) rho1,sigma(1),sigma(2)
c881   format(3f10.6)
      
      ymar(iwq)=(dt(nob+1,2)-pred(2))/sigma(2)  ! marginal 
      yjoint(iwq)=(dt(nob+1,1)-pred(1)-sigma(1)*ymar(iwq)*rho1)/rho11

      if(pred(2).lt.dt(nob+1,2)) iprob2=iprob2+1
      if(pred(1).lt.dt(nob+1,1)) iprob1=iprob1+1
      write(*,*) iprob1,iprob2
       
      
c****   product distribution function   ******
      xx1=ymar(iwq)
      xx2=yjoint(iwq)       
      write(17,6668) xx1	,xx2
6668  format(2f20.10)
      
      zmar=anordf(xx1) 
      zjoint=anordf(xx2)
       
      pcdf(iwq)=zmar*zjoint
      pcdf(iwq)=pcdf(iwq)-pcdf(iwq)*log(pcdf(iwq))
      
c****     ratio distribution function            
      cdfr(iwq)=zjoint/zmar
c      write(17,*) 'zjoint', zjoint,zmar,cdfr(iwq),iwq
      
      if (cdfr(iwq).lt.1.0) then
      cdfr(iwq)=0.5*cdfr(iwq)
      else
      cdfr(iwq)=1.0-(0.5/cdfr(iwq)) 
      endif     

99999 continue ! end big loop iwq=1,...,s 

      write(17,*)
      write(17,*) 'One-step ahead forecasting: '
      write(17,*) 'First set: product-smooth test'
      write(17,*) 'Second set: ratio '

      call smooth(pcdf,s,uu1,sprob)
      write(iout,2001) uu1(1),uu1(2),uu1(3),uu1(4),uu1(5)
2001  format(5f15.8)

      call smooth(cdfr,s,uu2,sprob)
      write(iout,2002) uu2(1),uu2(2),uu2(3),uu2(4),uu2(5)
2002  format(5f15.8)
c
      write(17,*)
      write(17,*) 'Kolmogorov-Smirnov test: 1st product, 2nd ratio'
      
      call ksonen(pcdf,s,uu3,sprob) ! KS-test product
      call ksonen(cdfr,s,uu4,sprob) !         ratio      
      write(iout,2003) uu3,uu4
2003  format(2f15.8)
c
       print*,'Output file: fort.17'
       stop
       end
cccccccccccccccccccccccccccccccccccccccccccccccccc
c     SUBROUTINES
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mtinv(a, da, kcomp,idim)
c**** 
      integer i,j,kcomp,idim,ii
      real*8 a(kcomp,kcomp), da(kcomp,kcomp)
      real*8 dpivot, pivot, deta, t

      deta = 1.0d0
c**** 
      if (idim .eq. 1) goto 600
      do 100 i = 1, idim
      do 100 j = 1, idim
  100 da(i,j) = a(i,j)
  120 do 500 i = 1, idim
      pivot = da(i,i)
c**** 
c**** DIVIDE PIVOT ROW BY PIVOT ELEMENT.
c**** 
      deta = deta * pivot
      da(i,i) = 1.0d0
      dpivot = pivot + 1.0d-25
      dpivot = da(i,i) / dpivot
      pivot = dpivot
      do 200 j = 1, idim
c**** 
c**** REDUCE NON-PIVOT ROWS
c****
  200 da(i,j) = da(i,j) * pivot
  210 do 500 ii = 1, idim
      if (ii .eq. i) goto 500
      t = da(ii,i)
      da(ii,i) = 0.0d0
      do 300 j = 1, idim
c****
  300 da(ii,j) = da(ii,j) - (da(i,j) * t)
  500 continue
      return 
  600 da(1,1) = deta / a(1,1)
      return 
      end
c****
c**** END OF 'MTINV'
c**** 
c*****************************************
       subroutine regr(y,x,maxn,md,ip,n,xpx,xpxinv,xpy,phi)
c
       real*8 y(*),x(maxn,md),xpx(md,md),xpy(*),phi(*)
       real*8 xpxinv(md,md), tem
c
       integer ip,n,i,j,it
c
       do 100 i=1, ip
        do 90 j = 1, i
         tem = 0.0d0
         do 80 it = 1, n
 80         tem=tem+x(it,i)*x(it,j)
         xpx(i,j) = tem
 90      xpx(j,i) = tem
 100     continue
c
        do 110 i=1, ip
         tem = 0.0d0
         do 105 it=1, n
 105       tem=tem+x(it,i)*y(it)
 110       xpy(i) = tem
c
        call mtinv(xpx,xpxinv,md,ip)
c
        do 150 i=1, ip
         tem = 0.0d0
         do 140 j = 1, ip
 140        tem=tem+xpxinv(i,j)*xpy(j)
 150        phi(i) = tem
c
        return
        end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c*****************************************
       subroutine detmtx(a,wk,mx,k,det)
c
c---- This program computes the determinant of a covariance matrix.
c     It assumes that the variance of each variable is not close to zero.
c     The program transforms the mtx into an upper triangular one.
c
       real*8 a(mx,mx), det, tem, wk(mx,mx), crit
       integer mx, k, i, j, km1
c--------- wk(.,.) is a working area.
c
       if(k.eq.1)then
        det = a(1,1)
        return
       endif
c
       km1 = k-1
       crit = 1.0d-15
c
       do 10 i=1, k
        do 10 j=1, k
 10        wk(i,j) = a(i,j)
c
       do 50 i=1, km1
        if(dabs(a(i,i)).le.crit)then 
         det = 0.0d0
         return
        endif
        tem = wk(i+1,i)/a(i,i)
        do 40 j=i+1, k
 40       wk(i+1,j) = wk(i+1,j)-wk(i,j)*tem
 50      continue
c
       det = 1.0d0
       do 60 i=1, k
 60       det = det*wk(i,i)
c
cc       print*,'det = ', det
c
       return
       end
ccccccccccccccccccccccccccccccccccccccccccccccccccc
C  
      REAL*8 FUNCTION RANDU(IX)  
C  
C UNIFORM RANDOM NUMBER GENERATOR.  
C REF: SCHRAGE, L. 1979 'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR'  
C      A. C. M. TRANSACTIONS ON MATHEMATICAL SOFTWARE V5 132-138.  
C
      implicit real*8 (a-h,o-z)   
      INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K  
      DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/  
      XHI=IX/B16  
      XALO=(IX-XHI*B16)*A  
      LEFTLO=XALO/B16  
      FHI=XHI*A+LEFTLO  
      K=FHI/B15  
      IX=(((XALO-LEFTLO*B16)-P)+(FHI-K*B15)*B16)+K  
      IF (IX.LT.0) IX=IX+P  
      RANDU=real(IX)*4.656612875d-10  
      RETURN  
      END  
      
cc******************************************************
      subroutine mtarsch(nob,dat,gthr)
cc
cc  This program is intended for nonlinear time series analysis.
cc  It searches for the threshold values (regime = 3, at most) of a
cc  multivariate TAR model. It uses Minimum AIC criterion.
cc
cc  Created on May 20, 1997 to analyze intraday data sets.
cc
cc  The program estimates a threshold-regression model: threshold variable
cc  and ranges of the thresholds must be given.
cc  input: dat(.,.) is the data matrix
cc  locy: the column number in dat(.,.) that gives the dependent var.
cc  locthr: column in which threshold variable is stored in dat(.,.)
cc  locx: columns where independent variables are located in dat(.,.).
cc  iorx: order of each independent variables
cc  ixlag: lags of each independent variables
cc  p: order of the dependent variable
cc  nx: number of exogeneous variables.
cc  (nst, nend): data span used in the program
cc         
cccccccccccccc Current setting
c    max-variables = 5
cc   max-obs = 10000
c
c       parameter(maxn=10000,md=50,mk=5)
       real*8 xpx(50,50),xpy(50,5),xpxinv(50,50) 
       real*8 phi(50,5),y(10000,5),dat(10000,6), thr(10)
       real*8 resi(10000,5), tem, chk, x(10000,50)
       real*8 aic, dlog, wk(50,50)
       real*8 thrva(500,2), bnd(2,2), aicmin, gthr(10)
c
       integer nob, p, ip, lagth, nreg, npts(10)
       integer locy(5), locthr, nx, ny, locx(5), nvar, i, j 
       integer iorx(5), ist, ia, ii, jj, it, i3, i4
       integer icnt, nst, nend, idx, i1, i2, npos(2)
       integer ixlag(10,5)
c-------- npos(.) stores the number of candidates for each threshold.
c
       character nfile*30
c
       iout = 15
       md=50
       maxn=10000
       mk=5
c
c---------- nob = # of observations; p = ar-order
c
c       write(6,1)
c    1  format(1x,'nob, nvar, n_y, n_x, delay, nregime, loc_thr: ',$)
c       read(5,*)nob, nvar, ny, nx, lagth, nreg, locthr

c       write(*,*) 'nob in mtarsch ', nob
       nvar=3  ! fixed values
       ny=2
       nx=1
       nreg=2
       locthr=3
       
c       if(nob .le. 0)stop
c       if(nob .gt. maxn)nob = maxn
c       if(nvar .le. 0)nvar = 1
c       if(nreg.le.1)stop
c       if(nreg.gt.3)then
c        print*,'too many regimes. Reset to 3'
c        nreg = 3
c       endif
       if((ny+nx).gt.nvar)then
        print*,'Number of variables fails to match'
        stop
       endif
c       if(nvax .gt. 5)then
c        print*,'too many variables! (the maximum is 5)'
c        stop
c       endif
       if(lagth.lt.0)lagth = 0
       if(locthr.le.0)locthr = 1
c
c       write(6,2)
c    2  format(1x,'Data_file Name: ',$)
c       read(5,3)nfile
c    3  format(a30)
c       open(unit=21,file=nfile,status='old')
c        write(*,*) 'nob and nvar', nob,nvar
c        do 4 i = 1, nob
c  4        write(*,*)(dat(i,j),j=1,nvar),nob
c       close(21)
c
c       write(6,401)
c 401   format(1x,'input data span (1 = yes): ',$)
c       read(5,*) ii
       ii=0
       if(ii.eq.1)then
        write(6,402)
        read(5,*) nst, nend
        if(nend.gt.nob) nend=nob
        if(nst.lt.1) nst = 1
        nob = nend-nst+1
        if(nst.eq.1)go to 415
        do 410 it=1, nob
         do 408 i=1, nvar
 408        dat(it,i) = dat(it-1+nst,i)
 410        continue
       endif
 402    format(1x,'data span (nst, nend): ',$)
 415    continue
c
c       write(6,5)
c 5     format(1x,'loc_y: ',$)
c       read(5,*) (locy(i),i=1,ny)
       locy(1)=1
       locy(2)=2
       locx(1)=3
c
c       if(nx.gt.0)then
c        write(6,7)
c        read(5,*)(locx(i),i=1,nx)
c       endif
c 7      format(1x,'loc_x: ',$)
c       write(6,8)
c 8     format(1x,'AR order: ',$)
c       read(5,*) p
       p=8  
       if(nx.gt.0)then
c        write(6,9)
c        read(5,*)(iorx(i),i=1,nx)

       iorx(1)=1  ! fixed  ****************
       
        do 44 i=1, nx
         do 42 j=1, iorx(i)
 42         ixlag(j,i) = j
 44         continue
c        write(6,56)
c        read(5,*) it
        it=0  !  fixed ***************
        
c        write(*,*) iorx(1),locy(1),locy(2)
        
        if(it.eq.1)then
         do 54 i=1, nx
c          write(6,57) i
c          read(5,*)(ixlag(j,i),j=1,iorx(i))
          do 53 j=1, iorx(i)
           if(ixlag(j,i).lt.0)then
              print*,'negative order not allowed'
              stop
             endif
 53        continue
 54       continue
        endif
       endif
c 9      format(1x,'#(lags) of the x-variables: ',$)
c 56     format(1x,'input detailed lags of exog. vars? (1=y): ',$)
c 57     format(1x,'lags of',i3,'-th exo. var (L to H): ',$)
c      
c--------- input thresholds
c       do 600 i=1, nreg-1
c       write(6,11) i
c       read(5,*)(bnd(i,j),j=1,2)
c 600   continue
c   11  format(1x,'input LOW. & UPP. bnd for threshold',i3,': ',$)
       bnd(1,1)=-.115  ! fixed bounds for thresholds *****
       bnd(1,2)=.145  
c       write(*,*) 'bnd ', bnd(1,1),bnd(1,2)
c
       thr(1) = -9.9d15
       i = nreg+1
       thr(i) = 9.9d15
c
c       write(6,610)
c 610   format(1x,'Input the numbers of pts for thresholds: ',$)
c       read(5,*)(npos(i),i=1,nreg-1)
       npos(1)=300   ! fixed  **************
       
       do 612 i=1, nreg-1
          if(npos(i).gt.500)then
           print*,'thr-npt exceeds 500, reset = 500'
           npos(i) = 500
          endif
 612      if(npos(i).le.0)npos(i) = 1
c
       do 650 i=1, nreg-1
        tem = (bnd(i,2)-bnd(i,1))/dfloat(npos(i)+1)
        do 645 j=1, npos(i)
 645       thrva(j,i) = bnd(i,1)+tem*dfloat(j)
c        write(*,*) 'tem = ', tem
 650       continue
c
      ist = p
      if(nx.gt.0)then
       do 17 i=1, nx
 17       if(ist.lt.ixlag(iorx(i),i))ist=ixlag(iorx(i),i)
      endif
      ist = max0(ist, lagth)+1
c
c------- including a constant term (ip = dimension of x-matrix)
       ip = 1+p*ny  
       if(nx.gt.0)then
        do 18 i=1, nx
 18        ip = ip + iorx(i)
       endif
       if(ip.gt.md)then
        print*,'dimension exceeds program specification'
        stop
       endif
c
      aicmin = 9.9d10
c--------------------------- prepare for the case of nreg = 2.
      if(nreg.eq.2)then
       npos(2) = 1
       thrva(1,2) = thr(3)
      endif
c----------------------------- 3 regimes are the maximum
      do 5000 i1 = 1, npos(1)
       thr(2) = thrva(i1,1)
       do 4500 i2 = 1, npos(2)
        thr(3) = thrva(i2,2)
        aic = 0.0d0
c
       do 3000  ia = 1, nreg
        icnt = 0
c
        do 100 it = ist, nob
         chk = dat(it-lagth,locthr)  
         if(chk.lt.thr(ia))go to 100
         if(chk.ge.thr(ia+1))go to 100
c------------- the observation is in the ia-th regime
         icnt = icnt+1
c         write(*,*) 'icnt = ', icnt
         do 20 i=1, ny
 20         y(icnt,i) = dat(it,locy(i))
         x(icnt,1) = 1.0d0
         idx = 1
         if(p.gt.0)then
          do 32 ii=1, ny 
           do 30 i=1, p
             idx=idx+1
 30          x(icnt,idx) = dat(it-i,locy(ii))
 32        continue
         endif
c
        if(nx.gt.0)then
         do 50 ii=1, nx
          if(iorx(ii).gt.0)then
           do 40 jj=1, iorx(ii)
              idx=idx+1
 40           x(icnt,idx)=dat(it-ixlag(jj,ii),locx(ii))
          endif
 50      continue
        endif
c
 100    continue
c
        npts(ia) = icnt
        if(icnt.le.ip)then
         print*,'insufficient data points in the', ia,'-th regime'
         go to 3000
        endif
c
       if(idx.ne.ip)then
        print*,'dimension error'
        stop
       endif
       call regrn(y,x,maxn,md,ny,ip,icnt,xpx,xpxinv,xpy,phi)
c
c----------residuals
       do 120 it=1, icnt
        do 115 i=1, ny 
         tem = 0.0d0
         do 110 ii=1, ip
 110      tem = tem + phi(ii,i)*x(it,ii)
       resi(it,i) = y(it,i)-tem
 115   continue
 120   continue
c----------------- Residual covariance matrix (MLE-type only for AIC).
       do 200 i=1, ny
        do 190 j=1, i
         tem = 0.0d0
         do 180 it=1, icnt
 180        tem=tem+resi(it,i)*resi(it,j)
        xpx(j,i) = tem/dfloat(icnt)
 190    xpx(i,j) = tem
 200    continue
c
       call detmtx(xpx,wk,md,ny,tem)
**
       tem = dlog(tem)*dfloat(icnt)
       aic = aic + tem + 2.0d0*dfloat(ip*ny)
c
c--------------- return for another regime
 3000  continue
c
c       write(iout,3004) aic
       if(aicmin .ge. aic)then
        aicmin = aic
        do 500 j=1, nreg-1
 500       gthr(j) = thr(j+1)
       i3 = i1
       i4 = i2
       endif
c
 4500  continue
 5000  continue
c
       write(6,3002)(gthr(i),i=1,nreg-1) ! fixed to one threshold
 3002  format(1x,'Thresholds: ',5f12.6)
       write(6,3003) aicmin
 3003  format(1x,'Overall minimum AIC',e12.5)
       write(6,3007) i3, i4
c 3004  format(1x,f14.6)
c       write(31,3005)(thrva(j,1),j=1,npos(1))
c       if(nreg.eq.3) write(32,3005)(thrva(j,2),j=1,npos(2))
c 3005  format(1x,f14.7)
 3007  format(1x,'index of thr-values: ',2i5)
c
c       print*,'AIC in fort.15, 1st thr in fort.31 & 2nd in fort.32'
c
       end
ccccccccccccccccccccccccccccccccccccccccccc
       subroutine regrn(y,x,maxn,md,ny,ip,n,xpx,xpxinv,xpy,phi)
c
       real*8 y(maxn,*),x(maxn,md),xpx(md,md),xpy(md,*),phi(md,*)
       real*8 xpxinv(md,md), tem
c
       integer ip,n,i,j,it,ii,ny
c
       do 100 i=1, ip
        do 90 j = 1, i
         tem = 0.0d0
         do 80 it = 1, n
 80         tem=tem+x(it,i)*x(it,j)
         xpx(i,j) = tem
 90      xpx(j,i) = tem
 100     continue
c
       do 120 ii=1, ny
        do 110 i=1, ip
         tem = 0.0d0
         do 105 it=1, n
 105       tem=tem+x(it,i)*y(it,ii)
 110       xpy(i,ii) = tem
 120       continue
c
        call mtinv(xpx,xpxinv,md,ip)
c
       do 160 ii=1, ny
        do 150 i=1, ip
         tem = 0.0d0
         do 140 j = 1, ip
 140        tem=tem+xpxinv(i,j)*xpy(j,ii)
 150        phi(i,ii) = tem
 160        continue
c
        return
        end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE smooth(dat,n,u,prob)
      INTEGER n,i1,i2
      real df
      real dat(n),prob(5),u(5)
!     computation first component
        u(1)=0
	u(2)=0
	u(3)=0
	u(4)=0
      ifail=0 
c      write(*,* ) 'n smooth = ', n
      do 10 i1=1,n
c	write(*,*) 'dat(i1)', dat(i1),i1
c	write(17,*) dat(i1),i1
       u(1)=u(1)+sqrt(12.0)*(dat(i1)-0.5)
       u(2)=u(2)+sqrt(5.0)*(6*((dat(i1)-0.5)**2)-0.5)
       u(3)=u(3)+sqrt(7.0)*(20*((dat(i1)-0.5)**3)-3*(dat(i1)-0.5))
       u(4)=u(4)+210*((dat(i1)-0.5)**4)-45*(dat(i1)-0.5)**2+(9.0/8.0)
c      write(17,*) u(1),u(2),u(3),u(4),i1
 10   continue  
      df=1.0
      do 12 i2=1,4
        u(i2)=(u(i2)**2)
c        write(17,*) ' u(i2) first ', u(i2)
        u(i2)=u(i2)/n
c        write(17,*) 'u ', u(i2),i2,n
 12   continue
      u(5)=0
      do 14 i1=1,4
         prob(i1)=1.0-chidf(u(i1),df)
c         write(7,*) 'prob ', prob(i1),df,i1
         u(5)=u(5)+u(i1)
 14   continue  
      df=4.0
      prob(5)=1.0-chidf(u(5),df)             
      return
      END 
c
      SUBROUTINE ksonen(datn,n,d,prob)
      INTEGER n
      REAL d,datn(n),prob
c      EXTERNAL func
CU    USES probks,sort
      INTEGER j
      REAL dt,en,ff,fn,fo,probks
      call sort(n,datn)
      en=n
      d=0.
      fo=0.
      do 11 j=1,n
        fn=j/en
        ff=datn(j)
        dt=max(abs(fo-ff),abs(fn-ff))
        if(dt.gt.d)d=dt
        fo=fn
11    continue
      en=sqrt(en)
      prob=probks((en+0.12+0.11/en)*d)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *$'.
c      
      FUNCTION probks(alam)
      real probks,alam,EPS1,EPS2
      PARAMETER (EPS1=0.001, EPS2=1.d-8)
      INTEGER j
      double precision a2,fac,term,termbf
      a2=-2.*alam**2
      fac=2.
      probks=0.
      termbf=0.
      do 11 j=1,100
        term=fac*dexp(a2*j**2)
        probks=probks+term
        if(abs(term).le.EPS1*termbf.or.abs(term).le.EPS2*probks)return
        fac=-fac
        termbf=abs(term)
11    continue
      probks=1.
      return
      END
!  (C) Copr. 1986-92 Numerical Recipes Software &*'>#)!D!"!*1.
             
      SUBROUTINE sort(n,arr)
      INTEGER n,M,NSTACK
      REAL arr(n)
      PARAMETER (M=7,NSTACK=50)
      INTEGER i,ir,j,jstack,k,l,istack(NSTACK)
      REAL a,temp
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.M)then
        do 12 j=l+1,ir
          a=arr(j)
          do 11 i=j-1,1,-1
            if(arr(i).le.a)goto 2
            arr(i+1)=arr(i)
11        continue
          i=0
2         arr(i+1)=a
12      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        temp=arr(k)
        arr(k)=arr(l+1)
        arr(l+1)=temp
        if(arr(l+1).gt.arr(ir))then
          temp=arr(l+1)
          arr(l+1)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l).gt.arr(ir))then
          temp=arr(l)
          arr(l)=arr(ir)
          arr(ir)=temp
        endif
        if(arr(l+1).gt.arr(l))then
          temp=arr(l+1)
          arr(l+1)=arr(l)
          arr(l)=temp
        endif
        i=l+1
        j=ir
        a=arr(l)
3       continue
          i=i+1
        if(arr(i).lt.a)goto 3
4       continue
          j=j-1
        if(arr(j).gt.a)goto 4
        if(j.lt.i)goto 5
        temp=arr(i)
        arr(i)=arr(j)
        arr(j)=temp
        goto 3
5       arr(l)=arr(j)
        arr(j)=a
        jstack=jstack+2
        if(jstack.gt.NSTACK)pause 'NSTACK too small in sort'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *$'.
   