The RATS Software Forum

Code: Select all

*
* @INOUETEST(options) x start end
*
* This procedure implements the test proposed in Inoue, A. (2001), Testing fo Distributional Change in Time Series, 
* Econometric Theory, 17, 156-187. Note that currently this only implements the Kolmogorov-Smirnov version of the test.
*
* Parameters:
*  x     - the time series to be tested
*  start - start of sample
*  end   - end of sample
*
* Options:
*  [PRINT]/NOPRINT
*  PLOT/[NOPLOT]
*  grid=vector containing the grid over which the supremum should be calculated
*  niter=number of iterations in bootstrap used to compute critical values [100]
*  block=window size used in the bootstrap 
*
* Revision Schedule:
*     12/2017 Written by Jonas Dovern, Heidelberg University
*
********************************************************************************************************************

procedure inouetest x start end
   type ser x
   type int start end
   *
   option switch  print 1
   option switch  plot 0
   option vec     grid
   option int     niter 100
   option int     block
   *
   local int  startl endl ngrid nobs tt bb mm gg
   local rec  stats1
   local ser  indicator z dif sumser meanser prodser
   local vec  lgrid T1star
   local real T1 sum1 sum2 sum3 Fhat
   local index temp
   *
   *** Check sample
   INQUIRE(series=x) startl>>start endl>>end
   comp nobs=endl-startl+1

   *** Determine grid
   if .not.%defined(grid) {
      statistics(noprint,fractiles) x startl endl
      comp lgrid=||%FRACT01,%FRACT05,%FRACT10,%FRACT25,%FRACT75,%FRACT90,%FRACT95,%FRACT99||
   }
   else {
      comp lgrid=grid
   }
   comp ngrid=%rows(lgrid)

   *** Obtain window size for bootstrap
   if .not.%defined(block) {
      comp lblock = fix(exp(0.134+0.499*log(nobs)))
   }
   else {
      comp lblock=block
   }

   *** Construct statistic
   dim stats1(nobs-1,%rows(lgrid))
   * Loop over all grid steps to find supremum
   do gg=1,ngrid
      set indicator startl endl = x<=lgrid(gg)
      statistics(noprint) indicator startl endl
      comp sum2 = %MEAN*%NOBS
      * Loop over all time periods to find supremum
      do tt=1,nobs-1
         statistics(noprint) indicator startl startl+tt-1
         comp sum1 = %MEAN*%NOBS
         comp stats1(tt,gg)=abs((1./sqrt(nobs))*(sum1-(tt*1./nobs)*sum2))
      end do bb
   end do tt
   * Compute test statistic, most likely break date, and the grid value that maximizes T1
   comp %T1=%maxvalue(stats1)
   comp temp=%%maxpos(stats1)
   comp %BREAKDATE=startl+temp(1)-1
   comp %MAXGRID=lgrid(temp(2))

   *** Obtain critical values by simulation
   dim T1star(niter)
   ewise stats1(i,j)=%na
   infobox(action=define,progress,lower=1,upper=niter) "Simulation ..."
   do bb=1,niter
      infobox(action=modify,current=bb)
      set z startl endl = %ran(sqrt(1./lblock))
      do gg=1,ngrid
         * Compute Fhat
         set indicator startl endl = x<=lgrid(gg)
         statistics(noprint) indicator startl endl
         comp Fhat=%MEAN

         * Compute series of difference between indicator function and Fhat
         set dif startl endl = %real(indicator(t))-Fhat

         * Compute sum over block length
         mvstats(means=meanser,width=lblock) dif startl+lblock-1 endl
         set sumser startl endl-lblock = meanser(t+lblock-1)*lblock

         * Compute product of random variable and block sum
         set prodser startl endl = z(t)*sumser(t)

         * Second sum does not depend on mm
         statistics(noprint) prodser startl endl-lblock+1
         comp sum2 = %MEAN*%NOBS

         do mm=1,nobs-lblock
            statistics(noprint) prodser startl startl+mm-1
            comp sum1 = %MEAN*%NOBS
            comp stats1(mm,gg) = abs((1./sqrt(nobs))*(sum1-(mm/(nobs-lblock+1.))*sum2))
         end do mm
      end do gg
      comp T1star(bb)=%maxvalue(%xsubmat(stats1,1,nobs-lblock,1,ngrid))
   end do bb
   infobox(action=remove)

   *** Compute critical values
   comp frac = %fractiles(T1star,||0.99,0.95,.9||)

   *** This prints test outcomes to the screen
   if print {
      disp; disp "Running Inoue (2001) test ..."
      disp; disp "Weighted KS statistic is T1=" %strval(%T1,"*.##")
      disp "It is obtained for " %datelabel(%BREAKDATE) "and a grid value of" %strval(%MAXGRID,"*.##")
      disp "Block size is" lblock
      disp; disp "Critical values based on" %string(niter) "iterations are:"
      disp " 1%:" frac(1)
      disp " 5%:" frac(2)
      disp "10%:" frac(3)
   }
   
   *** This plots the time series and puts a vertical line at the most likely break date
   if plot {
      set vline startl endl = t==%BREAKDATE
      graph(grid=vline) 1; # x startl endl
   }
end procedure