##################################################################
#simulated critical values for the limiting null distribution#
##############################################################
set.seed(100)
dd<-c(-.499,-(12:1)*0.04,0,(1:12)*0.04,.499)
numd<-length(dd)
result<-NULL
n<-5000
reptim<-30000
stat<-rep(0,reptim)
L0<-0.15
U0<-0.85
Ln0<-as.integer(n*L0)
Un0<-as.integer(n*U0)
sublength<-Un0-Ln0+1

for (h in 1:numd)
{
d1<-dd[h]
for (j in 1:reptim)
{
data<-as.vector(simARMA0(n,H=d1+1/2))
substat<-rep(0,sublength)
for (k in Ln0:Un0)
{
mean1<-mean(data[1:k])
mean2<-mean(data[(k+1):n])
inter1<-cumsum(data[1:k])-(1:k)*mean1
inter2<-cumsum(data[(k+1):n])-(1:(n-k))*mean2
substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n)
}
stat[j]<-max(substat)
print(j)
}

result<-rbind(result,as.vector(quantile(stat,c(0.9,0.95,0.99))))
print(h)
}


#################################################
#long memory, reptim=10000
#n=5000, simARMA0 from longmemo package in R
#d=-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499
#################################################
> result

h=1
         [,1]     [,2]     [,3]
[1,] 20.09091 21.99803 25.45926
  [,1]     [,2]     [,3]
[1,] 20.05159 21.92415 25.17966

          [,1]     [,2]     [,3]
 [1,] 20.23206 22.04757 25.57672
 [2,] 20.13137 22.14904 25.89856
 [3,] 20.90437 23.09339 26.81124
 [4,] 21.27919 23.71602 28.34925
 [5,] 21.60397 23.93517 29.22708
 [6,] 22.43455 25.07212 30.48307
 [7,] 23.04724 25.84554 31.76960
 [8,] 23.79224 26.91413 33.22689
 [9,] 24.31893 27.81315 34.73245
[10,] 25.31452 29.06364 36.63289
[11,] 26.07800 29.93718 38.13267
[12,] 27.21747 31.36806 40.58119
[13,] 28.12480 32.22254 41.98237
[14,] 28.70614 33.42419 43.32156
[15,] 29.93047 34.52267 45.57731
[16,] 31.21991 36.46689 48.11292
[17,] 32.51348 38.05722 50.15767
 [2,] 33.71634 39.24952 51.55109
 [3,] 34.82833 40.43071 53.81783
 [4,] 36.28219 42.84610 56.25739
 [5,] 37.18398 44.04682 59.46327
 [6,] 38.10746 45.33738 60.13217
 [7,] 39.84096 47.17937 62.52489
 [8,] 40.58492 48.16548 65.38060
 [9,] 41.86039 49.65368 68.31884
[10,] 43.85135 52.92012 71.68123
[11,] 43.76838 51.58495 69.85238

h=26 [2,] 43.50805 51.91513 68.68221

h=27
       [,1]     [,2]     [,3]
[1,] 44.7115 53.50515 71.58486

     [,1]     [,2]     [,3]
[1,] 43.91084 52.57802 72.04208

#############################
#n=5000, reptim=10000
#############################
20.05159 21.92415 25.17966
20.13137 22.14904 25.89856
20.90437 23.09339 26.81124
21.27919 23.71602 28.34925
21.60397 23.93517 29.22708
22.43455 25.07212 30.48307
23.04724 25.84554 31.76960
23.79224 26.91413 33.22689
24.31893 27.81315 34.73245
25.31452 29.06364 36.63289
26.07800 29.93718 38.13267
27.21747 31.36806 40.58119
28.12480 32.22254 41.98237
28.70614 33.42419 43.32156
29.93047 34.52267 45.57731
31.21991 36.46689 48.11292
32.51348 38.05722 50.15767
33.71634 39.24952 51.55109
34.82833 40.43071 53.81783
36.28219 42.84610 56.25739
37.18398 44.04682 59.46327
38.10746 45.33738 60.13217
39.84096 47.17937 62.52489
40.58492 48.16548 65.38060
41.86039 49.65368 68.31884
43.50805 51.91513 69.68221
44.71150 53.50515 71.58486
############################


############################
#local whittle estimation
############################


locwhittle<-function(d)
{
n<-length(ts)
xv<-2*pi*(1:m)/n
yv<-fft(ts[n:1])
yv<-yv[2:(m+1)]
yv<-(Mod(yv))^2/2/pi/n
obj1<-mean(xv^(2*d)*yv)
obj2<-mean(log(xv))*2*d
return(log(obj1)-obj2)
}
############################################33
#not used here
##############################################
Mmake.long<-function(d,ARts)
#ARts is H number of time series of length n
#[n,H] size, each column is a time series
{
L <- 10000
coef <- rep(0,L)
coef[1] <- 1
for (k in 1:(L-1))
  {
   coef[k+1] <- prod((d-1+(1:k))/(1:k))
}
n<-dim(ARts)[1]
H<-dim(ARts)[2]
result<-matrix(0,n,H)
for (k in 1:n)
{
result[k,1:H]<-matrix(coef[1:k],1,k)%*%ARts[k:1,1:H]
}
return(result[(n/2+1):n,1:H])
}

#######################################
#AR(1) model with Gaussian innovations
#######################################
MAR<-function(n,H,rr)
{
inter<-matrix(0,n+30,H)
epsilon<-matrix(rnorm((n+30)*H,0,1),(n+30),H)
for (j in 1:(n+29))
{
inter[j+1,1:H]<-epsilon[j+1,1:H]+rr*inter[j,1:H]
}
return(inter[31:(n+30),1:H])
}

#######################################
#MA(1) model with Gaussian innovations
#######################################
MMA<-function(n,H,rr)
{
inter<-matrix(0,n+30,H)
epsilon<-matrix(rnorm((n+31)*H,0,1),(n+31),H)
for (j in 1:(n+30))
{
inter[j,1:H]<-epsilon[j+1,1:H]+rr*epsilon[j,1:H]
}
return(inter[31:(n+30),1:H])
}
####################################3
#####################################3

#read the critical values (90%, 95%, 99%)

cv<-read.table("cv3.txt")
CV<-matrix(0,27,3)
for (j in 1:3)
{
CV[,j]<-cv[,j]
}

################################
#simulation setup (SIZE)
#####################
n<-100
reptim<-5000

###################
#long memory
###################
set.seed(100)
d<-0.4
data<-matrix(0,n,reptim)
for (j in 1:reptim)
{
data[,j]<-simARMA0(n,d+1/2)
}


#################
#AR(1)
set.seed(100)
data<-MAR(n,reptim,-0.5)
#################

####################
#MA(1)
set.seed(100)
data<-MMA(n,reptim,0.5)
#########################3


#################
#main program
#################
dn<-NULL
m<-as.integer(n^(2/3))

for (j in 1:reptim)
{
ts<-data[,j]
opt3 <- nlminb(start=c(0),obj=locwhittle,lower=-0.499,upper=0.499) 
dn<-c(dn,opt3$par)
}


critval<-matrix(0,2,reptim)
seq1<-c(-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499)

for (j in 1:reptim)
{
index1<-which.max(as.numeric(dn[j]<seq1+0.0000000001))
if (index1==1)
{
critval[1,j]<-CV[index1,1]
critval[2,j]<-CV[index1,2]
}
if (index1>1)
{
#90%
critval[1,j]<-((-dn[j]+seq1[index1])*CV[index1-1,1]+(dn[j]-seq1[index1-1])*CV[index1,1])/(seq1[index1]-seq1[index1-1])
#95%
critval[2,j]<-((-dn[j]+seq1[index1])*CV[index1-1,2]+(dn[j]-seq1[index1-1])*CV[index1,2])/(seq1[index1]-seq1[index1-1])
}
}

###############
stat<-rep(0,reptim)
L0<-0.15
U0<-0.85
Ln0<-as.integer(n*L0)
Un0<-as.integer(n*U0)
sublength<-Un0-Ln0+1
decision1<-rep(0,reptim)
decision2<-rep(0,reptim)

for (j in 1:reptim)
{
data1<-data[,j]
substat<-rep(0,sublength)
for (k in Ln0:Un0)
{
mean1<-mean(data1[1:k])
mean2<-mean(data1[(k+1):n])
inter1<-cumsum(data1[1:k])-(1:k)*mean1
inter2<-cumsum(data1[(k+1):n])-(1:(n-k))*mean2
substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n)
}
stat[j]<-max(substat)
#90%
decision1[j]<-as.numeric(stat[j]>critval[1,j])
#95%
decision2[j]<-as.numeric(stat[j]>critval[2,j])
print(j)
}
#return(c(sum(decision1)/reptim,sum(decision2)/reptim,dn[1:100]))
#}


round(sum(decision1)/reptim*100,2)

round(sum(decision2)/reptim*100,2)


#####################
#Power
#####################
n<-100
reptim<-5000

#############################
#series with a change point
##############################
unit1<-c(rep(0,.25*n),rep(1,.75*n))
unit2<-c(rep(0,.5*n),rep(1,.5*n))
unit3<-c(rep(0,.8*n),rep(1,.2*n))

extraMat1<-matrix(rep(unit1,reptim),n,reptim)
extraMat2<-matrix(rep(unit2,reptim),n,reptim)
extraMat3<-matrix(rep(unit3,reptim),n,reptim)

###################
#long memory
###################
set.seed(100)
d<-0
data<-matrix(0,n,reptim)
for (j in 1:reptim)
{
data[,j]<-simARMA0(n,d+1/2)
}

dataALT<-data+extraMat2*1


#################
#AR(1)
set.seed(100)
data<-MAR(n,reptim,-0.5)

dataALT<-data+extraMat3*0.5

#################

####################
#MA(1)
set.seed(100)
data<-MMA(n,reptim,0.5)
dataALT<-data+extraMat3*2
#########################3


##############################

#################
#main program
#################
dn<-NULL
m<-as.integer(n^(2/3))

for (j in 1:reptim)
{
ts<-dataALT[,j]
opt3 <- nlminb(start=c(0),obj=locwhittle,lower=-0.499,upper=0.499) 
dn<-c(dn,opt3$par)
}


critval<-matrix(0,2,reptim)
seq1<-c(-0.499, -(12:1)*0.04,0,(1:12)*0.04,.499)

for (j in 1:reptim)
{
index1<-which.max(as.numeric(dn[j]<seq1+0.0000000001))
if (index1==1)
{
critval[1,j]<-CV[index1,1]
critval[2,j]<-CV[index1,2]
}
if (index1>1)
{
#90%
critval[1,j]<-((-dn[j]+seq1[index1])*CV[index1-1,1]+(dn[j]-seq1[index1-1])*CV[index1,1])/(seq1[index1]-seq1[index1-1])
#95%
critval[2,j]<-((-dn[j]+seq1[index1])*CV[index1-1,2]+(dn[j]-seq1[index1-1])*CV[index1,2])/(seq1[index1]-seq1[index1-1])
}
}

###############
stat<-rep(0,reptim)
L0<-0.15
U0<-0.85
Ln0<-as.integer(n*L0)
Un0<-as.integer(n*U0)
sublength<-Un0-Ln0+1
decision1<-rep(0,reptim)
decision2<-rep(0,reptim)

for (j in 1:reptim)
{
data1<-dataALT[,j]
substat<-rep(0,sublength)
for (k in Ln0:Un0)
{
mean1<-mean(data1[1:k])
mean2<-mean(data1[(k+1):n])
inter1<-cumsum(data1[1:k])-(1:k)*mean1
inter2<-cumsum(data1[(k+1):n])-(1:(n-k))*mean2
substat[k-Ln0+1]<-sqrt(n)*abs(mean1-mean2)/(sqrt(sum(inter1^2)+sum(inter2^2))/n)
}
stat[j]<-max(substat)
#90%
decision1[j]<-as.numeric(stat[j]>critval[1,j])
#95%
decision2[j]<-as.numeric(stat[j]>critval[2,j])
#print(j)
}
#return(c(sum(decision1)/reptim,sum(decision2)/reptim,dn[1:100]))
#}


#sum(decision1)/reptim
#sum(decision2)/reptim



round(sum(decision1)/reptim*100,2)

round(sum(decision2)/reptim*100,2)