I would like to share with you a procedure that implements the test for a distributional change in a time series proposed in Inoue, A. (2001), Testing fo Distributional Change in Time Series, Econometric Theory, 17, 156-187.
Feel free to inform me about any bugs that you find.
Best,
Jonas
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