Nothing
#' 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 < 2) {warning("Frequency is < 2, 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)
if (sum(Augmentations)==0) {augmen = NULL
data <- data.frame(y, regular_integration = .Diff(.Lag(x,1),Freq), seasonal_integration=.Diff(.Lag(x,Freq),1))
} else {
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 <- stats::ts(utils::tail(stats::diffinv(y1, lag=Freq, xi=stats::rnorm(Freq,0,10)), Length), frequency=Freq, start=c(2001,1)) # seasonal diffinverting y1 to get a regular+seasonally intgr. series
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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.