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