R/ocsb.R

Defines functions ocsb

Documented in ocsb

#' OCSB test
#' 
#' Test for seasonal unit root roots in a time series.
#' @param x time series
#' @param freq frequency to be tested
#' @param method "OLS" or "ML"
#' @param augmentations non-seasonal and seasonal order of the augmentations
#' @param nrun number of runs in monte carlo simulation
#' @param seed seed for monte carlo simulated based generation of null distribution
#' @details The null hypothesis of the OCSB is that a series contains a seasonal unit root. This is tested by a Dickey-Fuller type regression. The test regression has often to be augmented by autocorrelational terms to ensure white noise of the error terms. 
#' @details If seasonal lags are included and method='OLS' the test regression is calculated by OLS, so only the seasonal lags are included. If instead of 'OLS' method='ML' a seasonal AR model is calculated, which implies that high-order non-seasonal lags will be indirectly included as well (see Box and Jenkins, 1970). For seasonal augmentations, ML is quite a bit slower than OLS. The run time can be speeded up by reducing the number of runs of the monte carlo simulation (e.g. nrun=100).
#' @details Under the null hypothesis the test statistic follows a non-standard distribution and thus needs to be simulated. The number of runs and the seed can be changed.
#' @references Box, G. and G. Jenkins (1970). Time Series Analysis: Forecasting and Control. San Francisco: Holden-Day.
#' @references Osborn D.R., Chui A.P.L., Smith J., and Birchenhall C.R. (1988). Seasonality and the order of integration for consumption, Oxford Bulletin of Economics and Statistics 50(4):361-377.
#' @author Daniel Ollech
#' @examples teststat <- ocsb(ts(rnorm(70, 10,10), frequency=7), nrun=200)
#' check_residuals(teststat)
#' @export


ocsb <- function(x, method="OLS", augmentations=c(3,0), freq=NA,nrun=1000, seed=123) {
  
if(method != "OLS" & method != "ML") {stop("Please, choose either OLS or ML as the method")}  
  if(method == "ML" & nrun > 100 & augmentations[2]>0) {cat("This will take some time. Please consider using method='OLS' or reducing nrun")}
  

if(is.na(freq)) {
  if (any(class(x)=="ts")) {freq <- stats::frequency(x)} else {
    if (any(class(x)=="xts")) {freq <- freq_xts(x)} else {
      stop("Do not know the frequency of the time series.")
    }}}
if(freq <= 1) {warning("Frequency is <= 1, this series does not seem to be seasonal. Results will almost surely be misguiding.")}
  
  if(method=="ML" & augmentations[2] > 0 & freq>249) {stop("If method = 'ML', seasonal lags cannot be included for a series with this many observations per year. Set method='OLS'.")} 
  
if (length(augmentations)==1) {augmentations = c(augmentations, 0)}
  
len <- length(x)
Start <- stats::start(x)

if (method=="ML") {  
  .teststat <- function(x, Augmentations=augmentations, Freq=freq) {
  mod <- tryCatch(suppressWarnings(stats::arima(.Diff(.Diff(x,1),Freq), order=c(Augmentations[1],0,0), seasonal=list(order=c(Augmentations[2],0,0), period=Freq), xreg=cbind(regular_integration =.Diff(.Lag(x,1),Freq),seasonal_integration=.Diff(.Lag(x,Freq),1)))) , error=function(e) e)
  if(inherits(mod, "error")) {mod <- NA;t_krit1 <- NA; stop("Error in the ARIMA model estimation, please try other version of OCSB test")} else {
    t_krit <- suppressWarnings(utils::tail(mod$coef,1)/utils::tail(sqrt(diag(mod$var.coef)),1))}
    return(list(t_krit=t_krit, mod=mod))
  }
  
  an <- .teststat(x)
  t_krit <- an$t_krit
  mod <- an$mod
  } 
  

  if (method=="OLS") {  
    
    .teststat <- function(x, Augmentations=augmentations, Freq=freq) {
    y <- .Diff(.Diff(x,1),Freq)
    augmen <- data.frame(lapply(1:max(1,Augmentations[1], Augmentations[2]*Freq), function(x) .Lag(y,k=x)))
    augmen <- augmen[, c(seq_len(Augmentations[1]), seq_len(Augmentations[2])*Freq)]
    nc_aug <- ifelse(sum(Augmentations)<=1, sum(Augmentations), ncol(augmen))
    if (nc_aug > 1) {colnames(augmen) <- paste0("augmentation", seq_len(nc_aug))}
    data <- data.frame(y, augmen, regular_integration = .Diff(.Lag(x,1),Freq), seasonal_integration=.Diff(.Lag(x,Freq),1))

    mod <- tryCatch(summary(stats::lm(y ~ ., data=data)), error=function(e) e)
    if(inherits(mod, "error")) {mod <- NA;t_krit1 <- NA; stop("Error in the OLS model estimation, please try other version of OCSB test")} else {
      t_krit <-  mod$coefficients[nrow(mod$coefficients),3]
    }
    return(list(t_krit=t_krit, mod=mod))
    }
    
    an <- .teststat(x)
    t_krit <- an$t_krit
    mod <- an$mod
  } 

.sim_int <- function(Length, Freq) {
  e <- stats::rnorm(Length*7, 0, 1) # Simulating error terms
  y1 <- stats::diffinv(e, xi=100)   # diffinverting error terms to get a regular integrated series
  y <- utils::tail(stats::diffinv(y1, lag=Freq, xi=stats::rnorm(Freq,0,10)), Length) # seasonal diffinverting y1 to get a regular+seasonally intgr. series
  y <- stats::ts(100+3*scale(y)*(stats::rbinom(1, 1,0.5)*2-1), frequency=Freq, start=c(2001,1))
  return(y)
}

.getPval <- function(x, Freq=freq, Seed=seed, Augmentations=augmentations, Len=len, sstart=Start) {
 #simulation of series
 set.seed(Seed + x)

 sim_series <- NA
 while(any(is.na(sim_series))) {
 sim_series <-  .sim_int(Len, Freq)} 
sim_series <- stats::ts(100 + scale(as.numeric(sim_series)), start=sstart, frequency=Freq)
    
 # get test statistic
 t_krit <- .teststat(sim_series, Augmentations=Augmentations)$t_krit
 return(t_krit)
}

null_distr <- sapply(seq_len(nrun), function(x) .getPval(x, Freq=freq, Seed=seed, Augmentations=augmentations, Len=len, sstart=Start)) 

p_value <- sum(t_krit > null_distr)/nrun

names(t_krit) <- ""
  
out <- list(stat=t_krit, Pval=p_value, test=paste0("OCSB with augmentations (regular, seasonal): (", augmentations[1], ", ", augmentations[2], ")"), model=mod, freq=freq)
class(out) <- c("seasinttests")
return(out)
}

Try the seastests package in your browser

Any scripts or data that you put into this service are public.

seastests documentation built on May 2, 2019, 2:37 a.m.