program aic_main
!***************************************************************************
! FORTRAN90 code (coded by Guodong Li)
! File: aic_new.f90
!
! Simulation experiment 2 in Biometrika paper Li and Li (2011).
!
! Phi =0 (size); Phi \neq 0 (power)
! Model: SETARMA(2;1,1,1,1) with delay d=1
!
! Reference:
! Li, G. and Li, W.K. (2011).
!   Testing a linear time series model against its
!   threshold extension. Biometrika, 98(1), 243-250.
!   DOI: 10.1093/biomet/asq074.
!
! IMSL (double precision) subroutines: 
!   DRLSE  (Fits a multiple linear regression using LS) 
!   DLINDS (Computes the inverse of a real symmetric positive definite matrix)
!   dEQTIL (compute empirical quantiles) 
!   drnnor (Generate pseudorandom numbers from a standard normal distribution 
!           using an inverse CDF method)
!   rnset  (Initialize a random seed for use in the IMSL randomnumber
!           generators)
!   drnun  (Generate pseudorandom numbers from a uniform (0,1) distribution)
!*****************************************************************************

   USE imsl
   IMPlICIT NONE
   integer,parameter::offset=200,n=200,replication=1200,Bootstrap=1000,pp=10
   real*8  innovation((n+offset)*replication),e(n+offset),x(n+offset)
   real*8  y(n)
   real*8  like0,like1,temp
   real*8  dist(Bootstrap*n),diss(n),test_star(bootstrap)!,fisher0(2,2)
   real*8  test_s(replication) !,test_s1(replication)
   integer NMISS,count,inde(replication)
   real*8  QPROP(2), Q(2), XLO(2), XHI(2),p_value(replication)
   integer i,iseed,time,repl,distribution,error_index,j,ii,jj
   real*8  psi0,phi0
   real*8  yy(n-pp),xx(n-pp,pp),b(pp),sst,sse,aic,aic1,r,b1(2*pp)
   real*8  xx1(n-pp,2*pp),hat_r
   integer p_aic,hat_n,bb
   real*8  test_stat,fisher(2*pp,2*pp),fisher0(2*pp,2*pp),score1(2*pp),temp1,q1(2),temp2

   iseed=time()
   CALL rnset(iseed)
   call drnnor((n+offset)*replication,innovation)

   distribution=1 ! 0 - uniform; 1 - two points distribution; 2 - normal

!  inde=0

do repl=1,replication
!******* start  replication ********

   !*** 1. Setting parameters
   psi0=0.1d0; phi0=0.0d0
   error_index=0

   !*** 2. Generate the sample
   e(1)=innovation((repl-1)*(n+offset)+1); x(1)=e(1); 
   e(2)=innovation((repl-1)*(n+offset)+2)*sqrt(0.5d0+0.5d0*e(1)*e(1))
   x(2)=e(2)
   
   do i=3,n+offset
      e(i)=innovation((repl-1)*(n+offset)+i)*sqrt(0.5d0+0.5d0*e(i-1)*e(i-1))
      x(i)=psi0*x(i-1)+0.1d0*x(i-2)+e(i)
      if(x(i-1).le.0.0d0) x(i)=x(i)+phi0*(x(i-1)+x(i-2))
   end do
   
   do i=1,n
      y(i)=x(i+offset)
   end do

   qprop=(/0.2d0,0.8d0/)
   CALL dEQTIL (n, y, 2, QPROP, Q, XLO, XHI, NMISS) ! percentiles

   !*** 3. Fit a AR(0) to AR(pp) models
   !***    Using AIC to select the order, denoted by p_aic

   do i=pp+1,n
      yy(i-pp)=y(i)
      do j=1,pp
         xx(i-pp,j)=y(i-j)
      end do
   end do

   aic=1.0d10
  
   do i=1,10
      call mlr(yy,n-pp,xx,pp,i,b,sst,sse)
      aic1=dble(n-pp)*log(sse/dble(n-pp))+2.0d0*dble(i)
      if(aic1.lt.aic) then
         p_aic=i
	 like0=sse
	 aic=aic1
      end if
   end do

   !*** 4. Fit a threshold AR(p_aic) model

   do i=pp+1,n
      do j=1,pp
         xx1(i-pp,2*j)=xx(i-pp,j)
      end do
   end do
   like1=like0

   do ii=1,n
      if((y(ii).lt.Q(2)).and.(y(ii).gt.Q(1))) then  ! start "if" loop
         r=y(ii)
        
         do i=pp+1,n
            do j=1,pp
               temp=0.0d0
               if(y(i-1).lt.r) temp=y(i-j)
	       xx1(i-pp,2*j-1)=temp
            end do
         end do
   
        call mlr(yy,n-pp,xx1,2*pp,2*p_aic,b1,sst,sse)

        if(sse.lt.like1) then
           like1=sse
	   hat_r=r
!	   hat_n=ii
        end if
     
     end if  ! end of "if" loop
   end do    ! ii=1,n

   !*** 5. The test statistic
   test_stat=(like0-like1)/(like0/dble(n-pp))

   !*** 6. Using bootstrap to approximate the critical values
      call mlr(yy,n-pp,xx,pp,p_aic,b,sst,sse)

    do i=pp+1,n
       temp=0.0d0;
       do j=1,p_aic
	  temp=temp+b(j)*y(i-j)
       end do
       e(i)=y(i)-temp
    end do

!   write(*,*)temp2

    qprop=(/0.05d0,0.95d0/)
    CALL dEQTIL (n, y, 2, QPROP, Q1, XLO, XHI, NMISS) ! percentiles
    
     if((Q1(1).gt.Q(1)).or.(Q1(2).lt.Q(2))) then
        error_index=1
        go to 555
     end if

     i=time()
     call rnset(i)

!***** choice of permutation disturb: uniform distribution
   if(distribution.eq.0) then
      call drnun(Bootstrap*n,dist)
      dist=dist-0.5
      dist=2.0*sqrt(3.0)*dist
   end if

!***** choice of the disturb : two points distribution
   if(distribution.eq.1) then
      call drnun(Bootstrap*n,dist)
      dist=dist-0.5
      do i=1,Bootstrap*n
         if(dist(i).gt.0.0d0) then
            dist(i)=1.0d0
         else
            dist(i)=-1.0d0
         end if
      end do
   end if

!***** choice of the disturb : standard normal distribution
   if(distribution.eq.2) then
      call drnnor(Bootstrap*n,dist)
   end if

   test_star=-10.0d0;

   do ii=1,n
      if((y(ii).lt.Q(2)).and.(y(ii).gt.Q(1))) then  ! start "if" loop
        r=y(ii)
	fisher=0.0d0
	
	do i=pp+1,n
	   b1=0.0d0
	   do j=1,p_aic
              b1(j)=y(i-j)
	   end do
	   do j=p_aic+1,2*p_aic
	      if(y(i-1).lt.r) b1(j)=y(i-j+p_aic)
	   end do

	   do j=1,2*p_aic
	      do jj=1,2*p_aic
	         fisher(j,jj)=fisher(j,jj)+b1(j)*b1(jj)
	      end do
	   end do
	end do    ! end ii=pp+1,n

        call sub1(2*pp,fisher,p_aic,2*p_aic,fisher0)

	do BB=1,bootstrap ! start bootstrap loop
           do i=pp+1,n
              diss(i)=dist((BB-1)*n+i)
           end do	

	   score1=0.0d0
	 
	   do i=pp+1,n
	      do j=1,p_aic 
	         score1(j)=score1(j)+e(i)*y(i-j)*diss(i)
	      end do
	      do j=p_aic+1,2*p_aic
	         if(y(i-1).lt.r) score1(j)=score1(j)+e(i)*y(i-j+p_aic)*diss(i)
	      end do
	   end do

	   temp=0.0d0
	 
	   do i=1,2*p_aic
	      do j=1,2*p_aic
	         temp=temp+score1(i)*fisher0(i,j)*score1(j)
	      end do
	   end do

           temp=temp/(like1/dble(n-pp))
	   if(test_star(BB).lt.temp) test_star(BB)=temp
	end do  ! end bootstrap loop
     end if ! end "if" loop
  end do    ! end ii=1,n

  count=0
  
  do i=1,bootstrap
     if(test_stat.lt.test_star(i)) count=count+1
  end do
  
555  p_value(repl)=dble(count)/dble(bootstrap)
  inde(repl)=error_index
  write(*,*) repl,test_stat,p_value(repl),error_index
  test_s(repl)=test_stat

!******* the end of replication ********
end do

open(10,file='res.txt')
do i=1,replication
   write(10,*) inde(i),p_value(i),test_s(i)
end do

close(10)

end program aic_main

subroutine mlr(y,n,x,p_max,p,bb,sst,sse)
  USE imsl
  IMPlICIT NONE
  integer n,p,p_max,intercept
  integer i,j
  real*8 y(n),x(n,p_max),xx(n,p)
  real*8 b(p),bb(p_max),sst,sse

  do i=1,n
     do j=1,p
        xx(i,j)=x(i,j)
     end do
  end do

  CALL dRLSE (n, y, p, xx, n, 0, B, SST, SSE)
  do i=1,p
     bb(i)=b(i)
  end do

  return
end subroutine mlr

subroutine sub1(n,a,n1,n2,ainv) !***  n2>n1
  USE imsl
  IMPlICIT NONE
  integer n,i,j,n1,n2
  real*8 a(n,n),a1(n1,n1),a2(n2,n2),ainv1(n1,n1),ainv2(n2,n2),ainv(n,n)
  
  do i=1,n1
     do j=1,n1
        a1(i,j)=a(i,j)
     end do
  end do

  do i=1,n2
     do j=1,n2
        a2(i,j)=a(i,j)
     end do
  end do

  CALL DLINDS (n1, A1, n1, AINV1, n1)
  CALL DLINDS (n2, A2, n2, AINV2, n2)
  
  do i=1,n1
     do j=1,n1
        ainv2(i,j)=ainv2(i,j)-ainv1(i,j)
     end do
  end do

  ainv=0.0d0
  do i=1,n2
     do j=1,n2
        ainv(i,j)=ainv2(i,j)
     end do
  end do

  return
end subroutine sub1
!  aic=dble(n-p_max)*log(sse/dble(n-p_max))+2.0d0*dble(p)+dble(n-p_max)*(1.0d0+log(6.283d0))
!  return