R/gof_stat.R

In GenOU: Sequential Change-Point Tests for Generalized Ornstein-Uhlenbeck Processes

#'@title Function to estimate quantiles for a goodness-of-fit test for generalized Ornstein-Uhlenbeck process
#'
#'@description Function to calculate the quantiles of Cramer-von Mise and Kolmogorov-Smirnov statistics.
#'
#'@param X     observations
#'@param T1    last time of observation
#'@param N     number of observations on from on interval (0,T1]
#'@param p     number of cosine coefficients >=1
#'@param q     number of sine coefficients >=0
#'@return \item{out}{List of statistics (cvm and ks), estimated parameters, and pseudo-observations}
#'@examples
#' T1=20
#' N=500
#' data(X)
#' out = gof_stat(X,T1,N,2,0)
#'@export
gof_stat= function(X,T1,N,p,q)
{

  DeltaN = T1/N # time increment Delta t

  n1=length(X)
  n=n1-1
  Y =matrix( diff(X),ncol=1)/ sqrt(DeltaN)
  tt = (c(1:n)-1)* DeltaN


  phi=matrix(1,ncol=1,nrow=n)
  if(p>1)
  {
    for(k in 1:(p-1))
    {
      phi=cbind(phi,fcos(tt*k))
    }
  }
  psi=NULL
  if(q>0)
  {
    for(k in 1:q)
    {
      psi=cbind(psi,fsin(tt*k))
    }
  }
  base=cbind(phi,psi)

  Z = cbind(base,-X[-n1])*sqrt(DeltaN)

  YN = Y[1:N]
  ZN = Z[1:N,]
  S= t(ZN) %*% ZN /T1
  sigmaHat = sqrt(sum(YN^2))/sqrt(N)
  thetaN = coef(lm(YN ~0+ZN))
  res = (YN-ZN%*%thetaN)/sigmaHat

  u=sort(pnorm(res))
  tt0=c(1:N)/N
  tt = tt0-0.5/N
  cvm = 1/(12*N)+ sum((u-tt)^2)
  ks=sqrt(N)*max(abs(u-tt0),abs(u-tt0+1/N))
  out=list(cvm=cvm,ks=ks,thetaN=thetaN, sigmaHat=sigmaHat, S=S, u=u)
  out
}