swh:1:snp:3a44eb759780145deea094ac2a25c5049546a085
Tip revision: c52ce072e457ea91212cea8cbcace0daa866910b authored by Han Lin Shang on 07 October 2018, 07:10:02 UTC
version 5.1
version 5.1
Tip revision: c52ce07
T_stationary.R
T_stationary <- function(sample, L = 49, J = 500, MC_rep=1000, cumulative_var = .90, Ker1 = FALSE, Ker2 = TRUE, h = ncol(sample)^.5, pivotal=FALSE)
{
xrefine = N = ncol(sample)
trefine = nrow(sample)
if(Ker1)
{
K=function(x)
{
output = min(1, max(1.1-abs(x),0))
return(output)
}
}
if(Ker2)
{
K=function(x)
{
output = min(1, max(2-2*abs(x),0))
return(output)
}
}
basis = create.fourier.basis(c(0,1),L)
ld = length(cumulative_var)
h1 = h
X1_bar = rowMeans(sample)
mean_subtracted_X1 = sample - X1_bar
gamma_hat = list()
for(i in 0:(N-1))
{
temp_matrix = matrix(rep(0, trefine^2), trefine,trefine)
for(j in (i+1):N)
{
temp_matrix = temp_matrix + (mean_subtracted_X1[,j] %*% t(mean_subtracted_X1[,j-i]))
}
gamma_hat[i+1]=list(temp_matrix/N)
}
cov_sample1 = gamma_hat[[1]]
for(index in 1:(N-1))
{
cov_sample1 = cov_sample1 + K(index/h1) *(gamma_hat[[index+1]] + t(gamma_hat[[index+1]]))
}
Z_matrix = cov_sample1
e1 = list()
for(index in 1:L)
{
e1[index] = list(as.matrix(eval.basis(evalarg =(1:trefine)/trefine, basisobj = basis,Lfdobj=0)[,index]))
}
eigenvalues1 = (eigen(Z_matrix)$values)/trefine
D = matrix(0,L,L)
for(k in 1:L)
{
for(ell in 1:L)
{
Integrand_matrix = Z_matrix * (e1[[k]] %*% t(e1[[ell]]))
D[k,ell] = 1/(trefine^2)*sum(Integrand_matrix)
}
}
eigenpairs = eigen(D)
eigenvectors = eigenpairs$vec
eigenvectors2 = eigen(Z_matrix)$vectors
eigenvalues = eigenpairs$val
evals = eigenvalues
if(pivotal)
{
d = c(1:ld)
switch = 0
stoper = 1
spot = 1
while(switch==0)
{
while((sum(eigenvalues[c(1:spot)])/sum(eigenvalues)) < cumulative_var[stoper])
{
spot = spot+1
}
d[stoper] = spot
stoper = stoper+1
if(stoper == (length(d)+1))
{
switch = 1
}
}
T_N0=1:ld
for(r in 1:ld)
{
ds=d[r]
inp.matrix=matrix(0,ds,N)
eig.v.norm=((trefine)^.5)*eigenvectors2
for(j in (1:ds))
{
for(k in (1:N))
{
inp.matrix[j,k]=t(sample[,k])%*%(eig.v.norm[,j])/trefine
}
}
T_Nsum=rep(0,ds)
for(j in (1:ds))
{
s.0=sum(inp.matrix[j,(1:xrefine)])
for(x in (1:xrefine))
{
T_Nsum[j]=T_Nsum[j]+(1/xrefine)*((1/N^.5)*(sum(inp.matrix[j,(1:x)])-(x/xrefine)*s.0))^2
}
}
T_N0[r]=sum(T_Nsum/eigenvalues[1:ds])
}
T = vector(, MC_rep)
T_array = matrix(0, length(d), MC_rep)
lambda = eigenvalues
for(dd in 1:length(d))
{
for(k in 1:MC_rep)
{
z=rnorm(d[dd]*J)
tot=0
for(n in c(1:d[dd]))
{
sum1 = 0
sum1 =sum((z[c(((n-1)*d[dd]+1):((n-1)*d[dd]+J))]/(pi*c(1:J)))^2)
tot = tot+sum1
}
T_array[dd,k] = T[k] = tot
}
}
p_values=(1:ld)
for(dd in 1:length(d))
{
p_values[dd] = round(1 - ecdf(T_array[dd,])(T_N0[dd]),4)
}
cat("\n")
cat("Pivotal test of stationarity for a functional time series\n")
cat("\n")
cat("null hypothesis: the series is stationary\n")
cat("\n")
cat(paste("p-value = ", p_values, sep=""),"\n")
cat(paste("N (number of functions) = ", N, sep=""),"\n")
cat(paste("number of MC replications = ", MC_rep, sep=""))
cat("\n")
}
if(pivotal==FALSE)
{
T = vector(, MC_rep)
T_array = (1: MC_rep)
lambda = eigenvalues
d = min(c(length(which(lambda>0)),15))
for(k in 1:MC_rep)
{
z=rnorm(d*J)
tot=0
for(n in c(1:d))
{
sum1 = 0
sum1 =sum((z[c(((n-1)*d+1):((n-1)*d+J))]/(pi*c(1:J)))^2)
tot = tot+lambda[n]*sum1
}
T_array[k] = T[k] = tot
}
int = sum(((1/sqrt(N)) *((sample[,1])-(1/xrefine)*rowSums(sample)))^2/(xrefine * trefine))
for(x in 2:xrefine)
{
int = int + sum(((1/sqrt(N)) * (rowSums(sample[,1:x]) -(x/xrefine) * rowSums(sample)))^2/(xrefine * trefine))
}
statT_N = int
p_values=(1:ld)
for(dd in 1:length(d))
{
p_values = round(1 - ecdf(T_array)(statT_N),4)
}
cat("\n")
cat("Monte Carlo test of stationarity of a functional time series\n")
cat("\n")
cat("null hypothesis: the series is stationary\n")
cat("\n")
cat(paste("p-value = ", p_values, sep=""),"\n")
cat(paste("N (number of functions) = ", N, sep=""),"\n")
cat(paste("number of MC replications = ", MC_rep, sep=""))
cat("\n")
}
}