c This is a fortran program using NAG library to generate a number of 

c alternative nonlinear SETAR models (of the form used in J.I.M.F.

c paper by Krager & Kugler) and compares the forecasts of an estimated SETAR 

c model with all estimated coefficients (including the threshold parameter, r, 

c the delay parameter, d) with upto an AR(4) model. These alternative 

c specifications are compared using the Diebold and Mariano test statistic

c on the FMSEs as well as by comparing the number of times the alternative

c models correctly put the variable in the correct regime.



      INTEGER n,i,j,k,p2,nrep,n1,nl,tnrep,edd

      PARAMETER (n=525,n1=500,tnrep=20000,nrep=1000,nl=4)

      DOUBLE PRECISION mse1(50,nl+1),mse2(50),yf(n),

     & var1(50,nl+1),bias1(50,nl+1),var2(50),bias2(50),

     & tet(n+8),cno1(3),cno2(3),cno3(3),ethr(2),

     & diff(50,nl+1),diff2(50,nl+1),cdiff(50,nl+1),tratio(50,nl+1),

     & ppar(50,nl+1),mmar(50,nl+1),nnar(50,nl+1),pnar(50,nl+1),

     & pmar(50,nl+1),npar(50,nl+1),nmar(50,nl+1),mpar(50,nl+1),

     & mnar(50,nl+1),

     & ppsar(50),mmsar(50),nnsar(50),pnsar(50),pmsar(50),npsar(50),

     & nmsar(50),mpsar(50),mnsar(50),yfar(n,nl+1),

     & dno1(3),dno2(3),dno3(3),svr1(3),svr2(3),svr3(3)

      DOUBLE PRECISION et(n),ratio1(50,nl+1),beta(12),jj1

      DOUBLE PRECISION g05ddf,syf2(n),cno(3),vr(3),dno(3),

     & b1(3),b2(4),b3(1),sb1(5),sb2(5),sb3(5),ab1(5),ab3(5),

     & ab2(5),cloc1(50,3),nloc1(50,3),cloc2(50,3),nloc2(50,3),

     & cloc3(50,3),nloc3(50,3),sig(4),thr(2),vb1(5),vb2(5),vb3(5),

     & ncorr(3),cncorr(3),f2

      EXTERNAL G05DDF,G05CBF

C m is the sample size and the FMSE are reported in temp1 and temp2

      open(1,file='dm_us13.dat')

      write(1,*) 'Estimating the Coeffs. of the unknown SETAR model'

      write(1,*) 'Lower Regime Exchange Rate Forecasts' 

c From Krager and Kugler (1993) (coeff(beta), sigma(sig), threshold (thr))

c Coefficients for the FF/$US

c      b1(1)=-0.0007

c      b2(1)=-0.0003

c      b2(2)=0.2642

c      b2(3)=0.0939

c      b2(4)=0.1148

c      b3(1)=0.0014

c      sig(1)=0.0157d0

c      sig(2)=0.0162d0

c      sig(3)=0.012d0

c      sig(4)=0.0189d0

c      thr(1)=-0.0029

c      thr(2)=0.0172

c Coefficients for DM/$US

      b1(1)=-0.0058

      b1(2)=-0.1736

      b1(3)=-0.1764

      b2(1)=-0.0003

      b2(2)=0.2881

      b2(3)=0.1222

      b2(4)=0.1470

      b3(1)=-0.001

      sig(1)=0.0171d0

      sig(2)=0.0181d0

      sig(3)=0.015d0

      sig(4)=0.0182d0

      thr(1)=-0.0096

      thr(2)=0.0165

c Coefficients LL/$US

c      b1(1)=0.0006

c      b2(1)=-0.0014

c      b2(2)=0.3620

c      b2(3)=0.1746

c      b2(4)=0.2101

c      b3(1)=0.0004

c      sig(1)=0.0158d0

c      sig(2)=0.0154d0

c      sig(3)=0.0126d0

c      sig(4)=0.0195d0

c      thr(1)=0.0041

c      thr(2)=0.0161

c Coefficients for Y/$US

c      b1(1)=-0.0030

c      b2(1)=-0.0014

c      b2(2)=1.3002

c      b2(3)=0.2918

c      b2(4)=0.2917

c      b3(1)=0.0010

c      sig(1)=0.0140d0

c      sig(2)=0.0174d0

c      sig(3)=0.0107d0

c      sig(4)=0.014d0

c      thr(1)=-0.0008

c      thr(2)=0.0036

c Coefficients for SF/$US

c      b1(1)=0.0008

c      b2(1)=-0.0011

c      b2(2)=0.2724

c      b2(3)=0.0501

c      b2(4)=0.1585

c      b3(1)=-0.0025

c      sig(1)=0.0202d0

c      sig(2)=0.0230d0

c      sig(3)=0.0162d0

c      sig(4)=0.0215d0

c      thr(1)=-0.0159

c      thr(2)=0.0208

c

      jj1=0.

c

      do 2 i=1,n-n1

      do 3 j=1,nl+1

      mse1(i,j)=0.0

      var1(i,j)=0.0

      bias1(i,j)=0.0

      diff(i,j)=0.0

      diff2(i,j)=0.0

      cdiff(i,j)=0.0

      tratio(i,j)=0.0

      ppar(i,j)=0.0

      nnar(i,j)=0.0

      mmar(i,j)=0.0

      npar(i,j)=0.0

      nmar(i,j)=0.0

      pnar(i,j)=0.0

      pmar(i,j)=0.0

      mpar(i,j)=0.0

      mnar(i,j)=0.0

   3  continue

      mse2(i)=0.0

      var2(i)=0.0

      bias2(i)=0.0

      ppsar(i)=0.0

      nnsar(i)=0.0

      mmsar(i)=0.0

      npsar(i)=0.0

      nmsar(i)=0.0

      pnsar(i)=0.0

      pmsar(i)=0.0

      mpsar(i)=0.0

      mnsar(i)=0.0

   2  continue

      do 4 j=1,3

      do 4 i=1,n-n1+1

      nloc1(i,j)=0.0

      nloc2(i,j)=0.0

      nloc3(i,j)=0.0

   4  continue

      do 181 i=1,3

      dno(i)=0.0

      vr(i)=0.0

      svr1(i)=0.0

      svr2(i)=0.0

      svr3(i)=0.0

      dno1(i)=0.0

      dno2(i)=0.0

 181  dno3(i)=0.0

      do 182 i=1,5

      ab2(i)=0.0

      ab1(i)=0.0

  182 ab3(i)=0.0

c

      call g05cbf(0)

      do 301 k=1,tnrep

c

c This generates the 500 observations within sample for one of the exchange 

c rate models (FF/$US,LL/$US,SF/$US,Y/$US)

c      tet(0)=g05ddf(0.0d0,sig(1))

c      tet(1)=g05ddf(0.0d0,sig(1))

c      tet(2)=G05DDF(0.0D0,sig(1))

c      DO 7 i=3,n+2

c      if (tet(i-1).le.thr(1))

c     & tet(i)=b1(1)+g05ddf(0.0d0,sig(2))

c      if (tet(i-1).gt.thr(2))

c     & tet(i)=b3(1)+g05ddf(0.0d0,sig(4))

c      if (tet(i-1).gt.thr(1).and.tet(i-1).le.thr(2)) then

c          f2=0.0

c          do 95 j=1,3

c   95     f2=f2+b2(j+1)*tet(i-j)

c          tet(i)=b2(1)+f2+g05ddf(0.0d0,sig(3))

c      end if

c    7 CONTINUE

c DM/$US rate from Krager & Kugler (1993)

      tet(0)=g05ddf(0.0d0,sig(1))

      tet(1)=g05ddf(0.0d0,sig(1))

      tet(2)=G05DDF(0.0d0,sig(1))

      do 7 i=3,n+2

      if (tet(i-1).le.thr(1)) then

         f2=0.0

         do 93 j=1,2

  93     f2=f2+b1(j+1)*tet(i-j)

         tet(i)=b1(1)+f2+g05ddf(0.0d0,sig(2))

      end if

      if (tet(i-1).gt.thr(2)) then

          tet(i)=b3(1)+g05ddf(0.0d0,sig(4))

      end if

      if (tet(i-1).gt.thr(1).and.tet(i-1).le.thr(2)) then

           f2=0.0

           do 95 j=1,3

   95      f2=f2+b2(j+1)*tet(i-j)

          tet(i)=b2(1)+f2+g05ddf(0.0d0,sig(3))

      end if

    7 continue

c

       do 103 i=1,n

  103  et(i)=tet(i+2)

c This selects whether we are in an lower regime

c      if (et(n1).gt.thr(1)) goto 301

c This selects whether we are in a middle regime

c      if (et(n1).le.thr(1).or.et(n1).gt.thr(2)) goto 301

c This selects whether we are in an upper regime

c      if (et(n1).le.thr(2)) goto 301



      jj1=jj1+1.

      if (jj1.ge.nrep+1.) goto 322

C This bit estimates from an AR(0) to an AR(nl) nl=4

      do 10 p2=0,nl

c

      call regress(p2,n1,et,beta)

c

c This forecasts the AR model 



      call forecast(p2,n1,et,beta,yf)



c This calculates the FMSE values for the AR models

      do 12 j=n1+1,n

      yfar(j,p2+1)=yf(j)

      mse1(j-n1,p2+1)=mse1(j-n1,p2+1)+((100.*(et(j)-yf(j)))**2)

     & /(1.*nrep)

      bias1(j-n1,p2+1)=bias1(j-n1,p2+1)+100.*(yf(j)-et(j))/(1.*nrep)

  12  continue

  10  continue



c

c This estimates the SETAR model



      call setar(n1,et,edd,ethr,syf2,

     & sb3,sb1,sb2,ncorr,cloc1,cloc2,cloc3)



c This forms the FMSEs as well as the bias and the variance for the SETAR 

c model. In addition, it calculates the proportion of times the actual and 

c forecast values fall into the same regime.

      do 14 i=n1+1,n

      bias2(i-n1)=bias2(i-n1)+100.*(syf2(i)-et(i))/(1.*nrep)

      mse2(i-n1)=mse2(i-n1)+((100.*(et(i)-syf2(i)))**2)/(1.*nrep)

      if (et(i).gt.thr(2).and.syf2(i).gt.thr(2)) 

     & ppsar(i-n1)=ppsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.syf2(i).le.thr(1)) 

     & nnsar(i-n1)=nnsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & syf2(i).gt.thr(1).and.syf2(i).le.thr(2))

     & mmsar(i-n1)=mmsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).gt.thr(2).and.syf2(i).le.thr(1)) 

     & pnsar(i-n1)=pnsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).gt.thr(2).and.syf2(i).gt.thr(1).and.

     & syf2(i).le.thr(2)) 

     & pmsar(i-n1)=pmsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.syf2(i).gt.thr(2)) 

     & npsar(i-n1)=npsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.syf2(i).gt.thr(1).and.

     & syf2(i).le.thr(2))

     & nmsar(i-n1)=nmsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & syf2(i).gt.thr(2))

     & mpsar(i-n1)=mpsar(i-n1)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & syf2(i).le.thr(1))

     & mnsar(i-n1)=mnsar(i-n1)+1.0/(1.0*nrep)

   14 continue





c This calculates the proportion of times the actual and forecast values fall 

c into the same regime for the AR models.

      do 615 j=1,nl+1

      do 616 i=n1+1,n

      if (et(i).gt.thr(2).and.yfar(i,j).gt.thr(2)) 

     & ppar(i-n1,j)=ppar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.yfar(i,j).le.thr(1)) 

     & nnar(i-n1,j)=nnar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & yfar(i,j).gt.thr(1).and.yfar(i,j).le.thr(2))

     & mmar(i-n1,j)=mmar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).gt.thr(2).and.yfar(i,j).le.thr(1)) 

     & pnar(i-n1,j)=pnar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).gt.thr(2).and.yfar(i,j).gt.thr(1).and.

     & yfar(i,j).le.thr(2)) 

     & pmar(i-n1,j)=pmar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.yfar(i,j).gt.thr(2)) 

     & npar(i-n1,j)=npar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).le.thr(1).and.yfar(i,j).gt.thr(1).and.

     & yfar(i,j).le.thr(2))

     & nmar(i-n1,j)=nmar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & yfar(i,j).gt.thr(2))

     & mpar(i-n1,j)=mpar(i-n1,j)+1.0/(1.0*nrep)

      if (et(i).gt.thr(1).and.et(i).le.thr(2).and.

     & yfar(i,j).le.thr(1))

     & mnar(i-n1,j)=mnar(i-n1,j)+1.0/(1.0*nrep)

c

c This forms the difference between the FMSE for the AR and SETAR models and

c then tests whether this difference is significant (Diebold and Mariano)



      diff(i-n1,j)=diff(i-n1,j)+((100.*et(i)-100.*syf2(i))**2-

     & (100.*et(i)-100.*yfar(i,j))**2)

      diff2(i-n1,j)=diff2(i-n1,j)+((100.*et(i)-100.*syf2(i))**2-

     & (100.*et(i)-100.*yfar(i,j))**2)**2

      do 617 ll=2,i-n1

      cdiff(i-n1,j)=cdiff(i-n1,j)+((100.*et(i-(ll-1))

     & -100.*syf2(i-(ll-1)))**2-

     & (100.*et(i-(ll-1))-100.*yfar(i-(ll-1),j))**2)*

     & ((100.*et(i)-100.*syf2(i))**2-

     & (100.*et(i)-100.*yfar(i,j))**2)

 617  continue

 616  continue

 615  continue





      do 16 j=1,3

      do 16 i=n1+2,n+1

      nloc1(i-n1,j)=nloc1(i-n1,j)+cloc1(i-n1,j)/(1.*nrep)

      nloc2(i-n1,j)=nloc2(i-n1,j)+cloc2(i-n1,j)/(1.*nrep)

      nloc3(i-n1,j)=nloc3(i-n1,j)+cloc3(i-n1,j)/(1.*nrep)

  16  continue

C This calculates the mean and variance for the number of observations in the

c 3 regimes.



      call summ(n1,edd,ethr,et,cno,cno1,cno2,cno3)



      do 22 i=1,3

      dno(i)=dno(i)+cno(i)/(1.*nrep)

      vr(i)=vr(i)+cno(i)*cno(i)/(1.*nrep)

      dno1(i)=dno1(i)+cno1(i)/(1.*nrep)

      dno2(i)=dno2(i)+cno2(i)/(1.*nrep)

      dno3(i)=dno3(i)+cno3(i)/(1.*nrep)

      svr1(i)=svr1(i)+(cno1(i)**2)/(1.*nrep)

      svr2(i)=svr2(i)+(cno2(i)**2)/(1.*nrep)

      svr3(i)=svr3(i)+(cno3(i)*cno3(i))/(1.*nrep)

  22  cncorr(i)=cncorr(i)+ncorr(i)/(1.*nrep)



C This calculates the mean and variance for the beta parameters in the SETAR 

c model



      do 24 i=1,5

      ab1(i)=ab1(i)+sb1(i)/(1.*nrep)

      vb1(i)=vb1(i)+(sb1(i)**2)/(1.*nrep)

      ab3(i)=ab3(i)+sb3(i)/(1.*nrep)

      vb3(i)=vb3(i)+(sb3(i)**2)/(1.*nrep)

      ab2(i)=ab2(i)+sb2(i)/(1.*nrep)

  24  vb2(i)=vb2(i)+(sb2(i)**2)/(1.*nrep)

 301  continue

 322  continue



      do 26 i=1,4

      vr(i)=vr(i)-dno(i)*dno(i)

      svr1(i)=svr1(i)-dno1(i)**2

      svr2(i)=svr2(i)-dno2(i)**2

 26   svr3(i)=svr3(i)-dno3(i)**2

      do 47 i=1,5

      vb1(i)=vb1(i)-ab1(i)*ab1(i)

      vb3(i)=vb3(i)-ab3(i)*ab3(i)

  47  vb2(i)=vb2(i)-ab2(i)*ab2(i)



c This calculates the ratio of the FMSEs of the AR to the SETAR models



      do 19 i=n1+1,n

      var2(i-n1)=mse2(i-n1)-(bias2(i-n1)**2)

      do 18 j=1,nl+1

      diff(i-n1,j)=diff(i-n1,j)/(1.0*nrep)

      cdiff(i-n1,j)=(diff2(i-n1,j)+2.*cdiff(i-n1,j))/(1.*nrep)

     & -diff(i-n1,j)**2

      var1(i-n1,j)=mse1(i-n1,j)-(bias1(i-n1,j)**2)

      ratio1(i-n1,j)=mse1(i-n1,j)/mse2(i-n1)

      tratio(i-n1,j)=diff(i-n1,j)/sqrt(cdiff(i-n1,j)/(1.0*nrep))

  18  continue

  19  continue



c This prints out the reslevant results

      write(1,*) 'FMSE of AR models relative to SETAR model'

      write(1,*) '   t      AR(0)     AR(1)     AR(2)     AR(3)

     & AR(4)'

      do 20 i=1,n-n1

      write(1,1011) 1.*i,(ratio1(i,j),j=1,5)

  20  continue

      write(1,*) 'T-ratio of AR(p) vs SETAR Model'

      do 204 i=1,n-n1

      write(1,1011) 1.*i,(tratio(i,j),j=1,5)

  204  continue

      write(1,*) 'Proportion of Times in A>0.0 & E>0.0'

      do 205 i=1,n-n1

      write(1,1019) 1.*i,(ppar(i,j),j=1,5),ppsar(i)

  205  continue

      write(1,*) 'Proportion of Times in A<0.0 & E<0.0'

      do 206 i=1,n-n1

      write(1,1019) 1.*i,(nnar(i,j),j=1,5),nnsar(i)

  206  continue

      write(1,*) 'Proportion of Times in A=M & E=M'

      do 207 i=1,n-n1

      write(1,1019) 1.*i,(mmar(i,j),j=1,5),mmsar(i)

  207  continue

      write(1,*) 'Proportion of Times in A>0.0 but E<0.0'

      do 208 i=1,n-n1

      write(1,1019) 1.*i,(pnar(i,j),j=1,5),pnsar(i)

  208  continue

      write(1,*) 'Proportion of Times in A>0.0 but E=M'

      do 209 i=1,n-n1

      write(1,1019) 1.*i,(pmar(i,j),j=1,5),pmsar(i)

  209  continue

      write(1,*) 'Proportion of Times in A<0.0 but E>0.0'

      do 210 i=1,n-n1

      write(1,1019) 1.*i,(npar(i,j),j=1,5),npsar(i)

  210  continue

      write(1,*) 'Proportion of Times in A<0.0 but E=M'

      do 211 i=1,n-n1

      write(1,1019) 1.*i,(nmar(i,j),j=1,5),nmsar(i)

  211  continue

      write(1,*) 'Proportion of Times in A=M but E<0.0'

      do 212 i=1,n-n1

      write(1,1019) 1.*i,(mnar(i,j),j=1,5),mnsar(i)

  212  continue

      write(1,*) 'Proportion of Times in A=M but E>0.0'

      do 213 i=1,n-n1

      write(1,1019) 1.*i,(mpar(i,j),j=1,5),mpsar(i)

  213  continue

      write(1,*) 'Bias of AR and SETAR models'

      do 23 i=1,n-n1

  23  write(1,1021) 1.*i,(bias1(i,j),j=1,5),bias2(i)

      write(1,*) 'Variance of AR and SETAR models'

      do 28 i=1,n-n1

  28  write(1,1021) 1.*i,(var1(i,j),j=1,5),var2(i)

      write(1,*) 'Proportion of times stay in the same regime'

      write(1,*) '       R1     R2      R3'

      write(1,*) '   t    2-periods               3-periods   

     &       4-periods'

      do 32 i=2,n-n1+1

  32  write(1,1013) 1.*(i-1),(nloc1(i,j),j=1,3),(nloc2(i,j),j=1,3),

     & (nloc3(i,j),j=1,3)

      write(1,*) 'Average Coefficient Estimates of SETAR model'

      write(1,1015) ab1(1),ab1(2),ab1(3),ab1(4),ab1(5),ab2(1),ab2(2),

     & ab2(3),ab2(4),ab2(5),ab3(1),ab3(2),ab3(3),ab3(4),ab3(5)

      write(1,*) 'SE of Coefficient Estimates of SETAR model'

      write(1,1015) sqrt(vb1(1)),sqrt(vb1(2)),sqrt(vb1(3)),

     & sqrt(vb1(4)),sqrt(vb1(5)),sqrt(vb2(1)),sqrt(vb2(2)),

     & sqrt(vb2(3)),sqrt(vb2(4)),sqrt(vb2(5)),sqrt(vb3(1)),

     & sqrt(vb3(2)),sqrt(vb3(3)),sqrt(vb3(4)),sqrt(vb3(5))

      write(1,*) 'Mean and Variance of % Occurance in 3 regimes'

      write(1,1012) dno(1),vr(1),dno(2),vr(2),dno(3),vr(3)

      write(1,*) 'Mean and Var. for % of time remain in regime'

      write(1,*) '    2-periods         3-periods          4-periods'

      write(1,1012) dno1(1),svr1(1),dno2(1),svr2(1),dno3(1),svr3(1)

      write(1,1012) dno1(2),svr1(2),dno2(2),svr2(2),dno3(2),svr3(2)

      write(1,1012) dno1(3),svr1(3),dno2(3),svr2(3),dno3(3),svr3(3)

      write(1,*) 'Proportion of times correct model selected'

      write(1,1014) cncorr(1),cncorr(2),cncorr(3)

 1011 format(f6.0,5f10.4)

 1012 format(6f9.3)

 1013 format(f3.0,9f8.3)

 1014 format(3f10.4)

 1015 format(6f10.4)

 1019 format(f6.0,6f10.4)

 1021 format(f6.0,6f12.6)

      STOP

      END

C

      subroutine regress(k1,n1,x,nb)

      INTEGER n1,i,ifail,k1,j,j1,nk

      parameter(n=525,nk=10)

      DOUBLE PRECISION y(n),nx(n,nk),xtx(nk+2,nk+2),x(n),

     & ixx(nk+2,nk+2),xy(nk+2),nb(nk+2),sigma2,tmp,res(n)

      EXTERNAL F01ADF,M01DAF

c This forms our x(i)=(1,y(i-1),y(i-2)..)' matrix and the y(i) vector

      do 2 i=1,nk+2

   2  nb(i)=0.0

      do 3 i=1,n1

      nx(i,1)=1.0

      y(i)=x(i)

   3  continue

      do 4 j=1,k1

      do 5 i=j+1,n1

      nx(i,(j+1))=x(i-j)

   5  continue

   4  continue

C This bit calculates the matrix x'x = xtx

      do 6 i=1,k1+1

      do 7 j=1,k1+1

      xtx(i,j)=0.0

      do 8 j1=k1+1,n1

      xtx(i,j)=xtx(i,j)+(nx(j1,i)*nx(j1,j))

   8  continue

   7  continue

   6  continue

      do 10 i=1,k1+1

      xy(i)=0.0

c This forms x'y vector

      do 10 j=k1+1,n1

  10  xy(i)=xy(i)+nx(j,i)*y(j)

      ifail=0

c This inverts the x'x matrix using the NAG library

      call f01adf(k1+1,xtx,nk+2,ifail)

      do 12 i=1,k1+1

      do 12 j=1,i

      ixx(i,j)=xtx(i+1,j)

  12  continue

C ixx is (X'X)^(-1)

      do 14 i=1,k1+1

      do 14 j=1,i-1

      ixx(j,i)=ixx(i,j)

   14 continue

c forms the OLS parameter estimates and writes this to nnb

      do 16 i=1,k1+1

      nb(i)=0.0

      do 16 j=1,k1+1

      nb(i)=nb(i)+ixx(i,j)*xy(j)

   16 continue

c      print *, nb(1)

      sigma2=0.0

c This calculates the residuals and sigma2

      do 18 i=k1+1,n1

      tmp=0.0

      do 20 j=1,k1+1

      tmp=tmp+nb(j)*nx(i,j)

  20  continue

      res(i)=y(i)-tmp

      sigma2=sigma2+(res(i)*res(i))/(1.*n1-2.*k1-1.)

  18  continue

      return

      end

C

      subroutine setreg(k1,n1,x,y1,sigma2,nbt)

      INTEGER n1,i,ifail,k1,j,j1,nk,n2

      parameter(n=525,nk=5)

      DOUBLE PRECISION y(n),nx(n,nk),xtx(nk+2,nk+2),x(n),

     & ixx(nk+2,nk+2),xy(nk+2),nbt(nk),sigma2,tmp,res(n),

     & y1(n),nnx(n,nk)

      EXTERNAL F01ADF,M01DAF

c This forms our x(i)=(1,y(i-1),y(i-2)..)' matrix and the y(i) vector

      do 23 j=1,nk

      do 2 i=1,n

   2  nnx(i,j)=0.0

  23  nbt(j)=0.0

      do 4 j=1,k1

      do 5 i=j+1,n1

      nnx(i,(j+1))=x(i-j)

   5  continue

   4  continue

      n2=0

      do 101 i=1,n1

      if (y1(i).ne.-999.99) then

          n2=n2+1

          y(n2)=y1(i)

          do 102 j=1,k1

 102      nx(n2,(j+1))=nnx(i,(j+1))

      endif

 101  continue

C This bit calculates the matrix x'x = xtx

      do 3 i=1,n2

   3  nx(i,1)=1.

      do 6 i=1,k1+1

      do 7 j=1,k1+1

      xtx(i,j)=0.0

      do 8 j1=1,n2

      xtx(i,j)=xtx(i,j)+(nx(j1,i)*nx(j1,j))

   8  continue

   7  continue

   6  continue

      do 10 i=1,k1+1

      xy(i)=0.0

c This forms x'y vector

      do 10 j=1,n2

  10  xy(i)=xy(i)+nx(j,i)*y(j)

      ifail=0

c This inverts the x'x matrix using the NAG library

      call f01adf(k1+1,xtx,nk+2,ifail)

      do 12 i=1,k1+1

      do 12 j=1,i

      ixx(i,j)=xtx(i+1,j)

  12  continue

C ixx is (X'X)^(-1)

      do 14 i=1,k1+1

      do 14 j=1,i-1

      ixx(j,i)=ixx(i,j)

   14 continue

c forms the OLS parameter estimates and writes this to nnb

      do 16 i=1,k1+1

      nbt(i)=0.0

      do 16 j=1,k1+1

      nbt(i)=nbt(i)+ixx(i,j)*xy(j)

   16 continue

c      print *, nbt(1),nbt(2)

      sigma2=0.0

c This calculates the residuals and sigma2

      do 18 i=1,n2

      tmp=0.0

      do 20 j=1,k1+1

      tmp=tmp+nbt(j)*nx(i,j)

  20  continue

      res(i)=y(i)-tmp

      sigma2=sigma2+(res(i)*res(i))/(1.*n2-2.*k1-1.)

  18  continue

      return

      end

c

      subroutine forecast(nl,n1,y,beta,yf)

      integer i,j,k,nl,n1,nob

      parameter(nob=525)

      DOUBLE PRECISION y(nob),beta(12),yf(nob),j1,j2

c Given our beta estimates for an AR model this forecasts the variable y(i)

      do 8 i=1,n1

      yf(i)=y(i)

  8   continue

      do 9 i=n1+1,nob

      yf(i)=0.0

  9   continue

      do 10 j=n1+1,nob

      j1=0.0

      j2=0.0

      do 11 k=1,nl

      j1=j1+beta(k+1)*yf(j-k)

  11  continue

      yf(j)=j1-j2+beta(1)

  10  continue

      return

      end

c

      subroutine setar(n2,dy,dd1,thres,

     & yf2,beta3,beta1,beta2,corr,loc1,loc2,loc3)

      integer n1,n2,i,j,r1,r2,nob,j1,j2,j3,jj,n1rep,dd,dd1,ncd,

     & cncd(3),ncncd(3),ifail

      parameter(nob=525,n1rep=500)

      double precision sr(nob),dy(nob),yf(nob),e1(nob),e2(nob),e3(nob),

     & beta1(5),beta2(5),beta3(5),be1(5),be2(5),be3(5),f2,

     & loc1(50,3),loc2(50,3),loc3(50,3),locl1(50),locl2(50),locl3(50),

     & yf2(nob),g05ddf,thres(2),thr(2),sigs(4),sig(4),

     & kb1(5),kb2(5),kb3(5),sigm1,sigm2,sigm3,aic,taic1,taic(3),laic(3),

     & corr(3),pr1(25)

      external g05ddf,m01caf



      open(2,file='dm_us11.dat')

      open(3,file='dm_us12.dat')

      n1=n2

      ifail=0

      do 3 i=1,n1

    3 sr(i)=dy(i)



c This sort the observations on sr from smallest to largest



      call m01caf(sr,1,n1,'ascending',ifail)



      taic1=100000.0



c The counter for the delay parameter

      do 2 dd=1,2



c the counter for the lower threshold value

      do 4 r1=1,71,2

      thr(1)=sr(n1*(r1+9)/100)



c The counter for the upper threshold value

      do 1 r2=r1+9,80,2

      thr(2)=sr(n1*(r2+10)/100)



      j1=0

      j2=0

      j3=0



c This assigns our variable to one of three sub-variables depending on y(i-1)

      do 6 i=1,n1

      e1(i)=-999.99

      e2(i)=-999.99

  6   e3(i)=-999.99

c dd is the delay parameter

          do 8 i=1+dd,n1

      if (dy(i-dd).le.thr(1)) then

          j1=j1+1

          e1(i)=dy(i)

      end if

      if (dy(i-dd).gt.thr(1).and.dy(i-dd).le.thr(2)) then

          j2=j2+1

          e2(i)=dy(i)

      end if

      if (dy(i-dd).gt.thr(2)) then

          j3=j3+1

          e3(i)=dy(i)

      end if

   8  continue



c This estimates our SETAR model as a conditional OLS with up to an 

c AR(5) in each regime and uses the AIC to choose the best model

c

      do 5 j=1,3

    5 corr(j)=0.0

      laic(1)=100000.00

      laic(2)=100000.00

      laic(3)=100000.00



      do 9 ncd=1,5

      call setreg(ncd-1,n1,dy,e1,sigm1,kb1)

      taic(1)=(1.*j1*log(sigm1))+2.*ncd

      if (taic(1).lt.laic(1)) then

          laic(1)=taic(1)

          do 71 j=1,ncd

   71     be1(j)=kb1(j)

          do 79 j=ncd+1,5

   79     be1(j)=0.0

          cncd(1)=ncd

          sigs(2)=sqrt(sigm1)

      end if

      call setreg(ncd-1,n1,dy,e2,sigm2,kb2)

      taic(2)=(1.*j2*log(sigm2))+2.*ncd

      if (taic(2).lt.laic(2)) then

          laic(2)=taic(2)

          do 72 j=1,ncd

   72     be2(j)=kb2(j)

          do 77 j=ncd+1,5

   77     be2(j)=0.0

          cncd(2)=ncd

          sigs(3)=sqrt(sigm2)

      end if

      call setreg(ncd-1,n1,dy,e3,sigm3,kb3)

      taic(3)=(1.*j3*log(sigm3))+2.*ncd

      if (taic(3).lt.laic(3)) then

          laic(3)=taic(3)

          do 73 j=1,ncd

   73     be3(j)=kb3(j)

          do 75 j=ncd+1,5

   75     be3(j)=0.0

          cncd(3)=ncd

          sigs(4)=sqrt(sigm3)

      end if

    9 continue



      aic=laic(1)+laic(2)+laic(3)

      if (aic.lt.taic1) then

          taic1=aic

          ncncd(1)=cncd(1)

          ncncd(2)=cncd(2)

          ncncd(3)=cncd(3)

          thres(1)=thr(1)

          thres(2)=thr(2)

          sig(2)=sigs(2)

          sig(3)=sigs(3)

          sig(4)=sigs(4)

          do 11 i=1,ncncd(1)

   11     beta1(i)=be1(i)

          do 81 i=ncncd(1)+1,5

   81     beta1(i)=0.0

          do 15 i=1,ncncd(2)

   15     beta2(i)=be2(i)

          do 83 i=ncncd(2)+1,5

   83     beta2(i)=0.0

          do 17 i=1,ncncd(3)

   17     beta3(i)=be3(i)

          do 85 i=ncncd(3)+1,5

   85     beta3(i)=0.0

          sig(1)=sqrt((j1*(sig(2)**2)+j2*(sig(3)**2)+j3*(sig(4)**2))

     &    /(j1+j2+j3))

          dd1=dd

      end if

    1 continue

c      print *, r1,r2,dd

    4 continue

    2 continue



      write(2,1012) 1.*dd1,thres(1),thres(2)

 1012 format(f6.1,2f12.8)



c NOTE: for DM/$US(p1=3,p2=4,p3=1) otherwise(p1=1,p2=4,p3=1)



      if (ncncd(1).eq.1.and.dd1.eq.1) corr(1)=1.0

      if (ncncd(2).eq.4.and.dd1.eq.1) corr(2)=1.0

      if (ncncd(3).eq.1.and.dd1.eq.1) corr(3)=1.0



      do 10 i=1,n1

  10  yf(i)=dy(i)

      do 101 i=1,nob

 101  yf2(i)=0.0

      do 102 i=n1+1,nob

 102  pr1(i-n1)=0.0

c This ensures we simulate over the regimes that can be observed.

      do 51 jj=1,n1rep

c This is to forecast upto nob-n1-1 observations ahead.

      loc1(0,1)=0.0

      loc1(1,1)=0.0

      loc1(0,2)=0.0

      loc1(1,2)=0.0

      loc1(0,3)=0.0

      loc1(1,3)=0.0

      do 12 i=n1+1,nob

      do 100 j=1,3

      loc1(i-n1+1,j)=0.0

      loc2(i-n1+1,j)=0.0

 100  loc3(i-n1+1,j)=0.0



c Using the parameter estimates above we forecast our variables conditional on

c yf(i-dd1). The loc(?) variable determines the number of times we 

c remain in the same regime for upto ? successive periods in the forecast

c period.



      if (yf(i-dd1).le.thres(1)) then

         f2=0.0

         do 14 j=1,ncncd(1)-1

  14     f2=f2+beta1(j+1)*yf(i-j)

         yf(i)=f2+beta1(1)+g05ddf(0.0d0,sig(2))

         if (yf(n1-2).le.thres(1)) loc1(0,1)=1.

         if (yf(n1-1).le.thres(1)) loc1(1,1)=1.

         if (yf(i).le.thres(1)) loc1(i-n1+1,1)=1.

      end if

      if (yf(i-dd1).gt.thres(1).and.yf(i-dd1).le.thres(2)) then

         f2=0.0

         do 16 j=1,ncncd(2)-1

  16     f2=f2+beta2(j+1)*yf(i-j)

         yf(i)=f2+beta2(1)+g05ddf(0.0d0,sig(3))

         if (yf(n1-2).gt.thres(1).and.yf(n1-2).le.thres(2)) 

     &    loc1(0,2)=1.

         if (yf(n1-1).gt.thres(1).and.yf(n1-1).le.thres(2))

     &    loc1(1,2)=1.

         if (yf(i).gt.thres(1).and.yf(i).le.thres(2)) 

     &    loc1(i-n1+1,2)=1.

      end if

      if (yf(i-dd1).gt.thres(2)) then

         f2=0.0

         do 20 j=1,ncncd(3)-1

  20     f2=f2+beta3(j+1)*yf(i-j)

         yf(i)=f2+beta3(1)+g05ddf(0.0d0,sig(4))

         if (yf(n1-2).gt.thres(2)) loc1(0,3)=1.

         if (yf(n1-1).gt.thres(2)) loc1(1,3)=1.

         if (yf(i).gt.thres(2)) loc1(i-n1+1,3)=1.

      end if

c

      loc2(i-n1+1,1)=loc1(i-n1+1,1)*loc1(i-n1,1)

      loc3(i-n1+1,1)=loc2(i-n1+1,1)*loc1(i-n1-1,1)

      loc2(i-n1+1,2)=loc1(i-n1+1,2)*loc1(i-n1,2)

      loc3(i-n1+1,2)=loc2(i-n1+1,2)*loc1(i-n1-1,2)

      loc2(i-n1+1,3)=loc1(i-n1+1,3)*loc1(i-n1,3)

      loc3(i-n1+1,3)=loc2(i-n1+1,3)*loc1(i-n1-1,3)

      locl1(i-n1+1)=loc1(i-n1+1,1)+loc1(i-n1+1,2)+loc1(i-n1+1,3)

      locl2(i-n1+1)=loc2(i-n1+1,1)+loc2(i-n1+1,2)+loc2(i-n1+1,3)

      locl3(i-n1+1)=loc3(i-n1+1,1)+loc3(i-n1+1,2)+loc3(i-n1+1,3)

C This section calculates the probability of observing a forecasting less

c than or equal to the actual value

      if (yf(i).lt.dy(i)) pr1(i-n1)=pr1(i-n1)+1./(1.*n1rep)

      yf2(i)=yf2(i)+yf(i)/(1.*n1rep)

  12  continue

  51  continue

      write(3,1111) pr1(1),pr1(2),pr1(5)

 1111 format(3f8.3)

      return

      end

c

      subroutine summ(n1,dd,thres,dy,co,co1,co2,co3)

      integer n,i,n1,j1,j2,j3,dd

      parameter (n=525)

      double precision dy(n),co(3),co1(3),co2(3),co3(3),thres(2)

      j1=0.0

      j2=0.0

      j3=0.0

      do 2 i=1,3

      co(i)=0.0

      co1(i)=0.0

      co2(i)=0.0

   2  co3(i)=0.0

c This assigns our variable to one of three sub-variables depending on y(i-dd)

c dd is the delay parameter and counts the number of times we are in each

c regime for 1, 2 and 3 successive periods over the entire sample.

      do 4 i=1+dd,n1

      if (dy(i-dd).le.thres(1)) then

          j1=j1+1

      end if

      if (dy(i-dd).gt.thres(1).and.dy(i-dd).le.thres(2)) then

          j2=j2+1

      end if

      if (dy(i-dd).gt.thres(2)) then

          j3=j3+1

      end if

   4  continue

      do 6 i=11+dd,n1

      if (dy(i-dd).le.thres(1)) then

        if (dy(i-dd-1).le.thres(1)) then

          co1(1)=co1(1)+1.0/((n1-10.-dd)/100.)

          if (dy(i-dd-2).le.thres(1)) then

            co2(1)=co2(1)+1.0/((n1-10.-dd)/100.)

            if (dy(i-dd-3).le.thres(1)) then

              co3(1)=co3(1)+1.0/((n1-10.-dd)/100.)

            end if

          end if

        end if

      end if

      if (dy(i-dd).gt.thres(1).and.dy(i-dd).le.thres(2)) then

        if (dy(i-dd-1).gt.thres(1).and.dy(i-dd-1).le.thres(2)) then

          co1(2)=co1(2)+1.0/((n1-10.-dd)/100.)

          if (dy(i-dd-2).gt.thres(1).and.dy(i-dd-2).le.thres(2)) then

            co2(2)=co2(2)+1.0/((n1-10.-dd)/100.)

            if (dy(i-dd-3).gt.thres(1).and.dy(i-dd-3).le.thres(2)) then

              co3(2)=co3(2)+1.0/((n1-10.-dd)/100.)

            end if

          end if

        end if

      end if

      if (dy(i-dd).gt.thres(2)) then

        if (dy(i-dd-1).gt.thres(2)) then

          co1(3)=co1(3)+1.0/((n1-10.-dd)/100.)

          if (dy(i-dd-2).gt.thres(2)) then

            co2(3)=co2(3)+1.0/((n1-10.-dd)/100.)

            if (dy(i-dd-3).gt.thres(2)) then

              co3(3)=co3(3)+1.0/((n1-10.-dd)/100.)

            end if

          end if

        end if

      end if

   6  continue

c These are the proportion of times in each regime

      co(1)=(j1*1.)/((n1-dd)/100.)

      co(2)=(j2*1.)/((n1-dd)/100.)

      co(3)=(j3*1.)/((n1-dd)/100.)



      return

      end

