[R] brillinger test

[vito muggeo]

Dear all,
below there a code to calculate the Brillinger statistic, testing for
non-monotonic trend in (time series) data.

The function (probably it is not optimized) works, but I could not check for
it. That is, I would be sure that it provides correct output. Does anyone
know any well-known data-set where I can check the output? Or is there
anyone being able to check for it?
Many thanks,
best,
vito


brillinger<-function(y,v=length(y)/10,L=length(y)/20,display=FALSE,iid=FALSE
){   #Brillinger test for non-monotonic trend. Brillinger (1989) Biometrika,
76:23-30.
#   require toeplitz() function in the ts package.
#y: the time series. It has to be a vector.
#v: parameter modelling the Moving Average-based trend estimate
#L: parameter modelling the variance estimate
#display: if TRUE the data are plotted with a nonparametric estimate
superimposed
#iid: are the observations iid?
#author: vito muggeo <vito.muggeo at giustizia.it>
    yy<-y
    y.lab<-deparse(substitute(y))
    n<-length(y)
    y<-y[(v+1):(n-v)]
    x<-1:length(y)
    c.t1<-sqrt((x-1)*(1-(x-1)/n))
    c.t2<-sqrt(x*(1-x/n))
    c.t<-c.t1-c.t2
    Y<-t(toeplitz(yy))
    Y[row(Y)<col(Y)]<-NA
    mu<-rowMeans(Y[,1:(2*v+1)],na.rm=FALSE)
    if(display) {plot(y,type="l",lty=3,xlab="Time",ylab=y.lab)
                lines(mu[!is.na(mu)])
                }
    e<-y-mu[!is.na(mu)]
        if(iid){z.iid<-sum(c.t*y)/(sd(e)*sqrt(sum(c.t^2)))
        p<-1-1*pnorm(abs(z.iid))
        out<-list("Brillinger Statistic"=z.iid,"p-value (one-sided)"=p)
        return(out)
        }
    csi<-apply(as.matrix(1:L),1,function(J){sum(e*exp(-2i*pi*J*(x-1)/n))})

a<-apply(as.matrix(1:L),1,function(J){sin(2*pi*J*(2*v+1)/(2*n))/((2*v+1)*sin
(2*pi*J/(2*n)))})
    den<-sum((1-a)^2)
    S0<-sum(Mod(csi)^2/n)/den
    z<-sum(c.t*y)/sqrt(S0*sum(c.t^2))
    p<-1-1*pnorm(abs(z))
    out<-list("Brillinger Statistic"=z,"p-value (one-sided)"=p)
    return(out)
    }

#Example
> n<-200
> y<-ifelse((1:n)>0.5*n,1,0)*2.5+rnorm(n) #non monotonic trend time serie
> library(ts)
> brillinger(y,display=TRUE,iid=TRUE)
$"Brillinger Statistic"
[1] 1.824676

$"p-value (one-sided)"
[1] 0.03402500