Gfrktyl

SHORT VERSION: Made a workflow function, want a better output with methods.

General Background

I have been working on a package for some time to eliminate the tedium of preliminary time series econometrics (unit root testing, cointegration, model building and lag/lead selection). I am aware of many alternatives out there but I am still happy with the direction I am going for now. Based off of Stock & Watson.

Function Workflow (simple version)

1. Takes a user-specified cointegrating relationship (written as a formula): $$Y_t = alpha_t + X_t$$ in R: Y ~ 1 + X
where the dependent and independent variables are all nonstationary and alpha is the (optional) constant.

2. Creates the formula
$$Y_t = alpha_t + X_t + sum_{i = -k}^k Delta X_{t-i}$$
in R: Y ~ 1 + X + L(diff(X),-k:k) using dynlm's handy L() capabilities (although I do not use D() for diff(). Where k is the maximum lag/lead value.

3. The rest of the function runs different k values and selects the best model using the BIC function (but any model selection function can be put in as an argument). It also computes HAC estimated errors using Newey West.

Code
Here is the function. I am very open to criticism/critiques/comments.

Requires packages:

library(dynlm)

library(lmtest)

library(sandwich)



buildDOLS <- function (coint_formula, data, fixedk = NULL, robusterrors = TRUE, selection = BIC){

  # checks

  stopifnot(is.ts(data)) # time series data

  stopifnot(is.null(fixedk)|is.numeric(fixedk)) # fixed k either is null or is numeric

  stopifnot(is.function(selection)) # selection method is a function (should work on a function)

  #  Formula creation

  ff <- coint_formula

  all_names <- dimnames(attr(terms(ff), "factors")) # X and Y variables

  y_names <- all_names[[1]][!(all_names[[1]] %in% all_names[[2]])]

  x_names <- all_names[[2]][all_names[[2]] %in% colnames(data)]

  # Dynamic Ordinary Least Squares formulation

  ff_LHS <- y_names

  ff_RHS <- paste(c(ifelse(attr(terms(ff), "intercept") == 1, "1", "-1"), # constant 

                    x_names, # input variables

                    paste0("L(diff(", x_names, "),-k:k)")), # sum of lead/lagged differences of x variables

                  collapse=" + ")

  ff_k <- paste(ff_LHS, "~", ff_RHS)

  # if k (the maximum number of lags/leads) was not fixed, use a default value

  k <- ifelse(is.null(fixedk), floor(dim(data)[1]^(1/3)/2), fixedk)

  # run the model. If k was fixed, this is the final model:

  DOLS_k <- dynlm(formula(ff_k), data = data) 

  # If k was not fixed, DOLS_k will be used to keep constant the start and end dates during model selection

  if(is.null(fixedk)){

    # Use any selection function that is indicated in the selection argument

    k_select <- sapply(1:k, function(k) match.fun(FUN = selection)(dynlm(formula(ff_k), 

                                                                         data = data, 

                                                                         start = start(DOLS_k), 

                                                                         end = end(DOLS_k))))

    # only re-estimate the model if k_select differs from k to be efficient

    if(k != which.min(k_select)){

      k <- which.min(k_select)

      DOLS_k <- dynlm(formula(ff_k), data = data)

    } 

    # save the selection matrix results inside the model

    DOLS_k$selection <- cbind(1:k, k_select, DOLS_k$df + length(DOLS_k$coeff), 

                              start(DOLS_k)[1], end(DOLS_k)[1])

    colnames(DOLS_k$selection) <- c("# of lags/leads (k)", deparse(substitute(selection)),"#Obs",

                                    "StartDate", "EndDate")

  }

  DOLS_k$k <- k # save the lag used inside the model

  # save the HAC estimated errors inside the model

  if(robusterrors) DOLS_k$HAC <- lmtest::coeftest(DOLS_k, vcov = sandwich::NeweyWest(DOLS_k, lag = k))

  # rewriting the call function to be a run on its own 

  DOLS_k$call <- as.call(c(quote(dynlm), 

                           formula = formula(gsub("-k:k", paste0("-",k,":",k), ff_k)), 

                           data = substitute(data)))                  

  class(DOLS_k) <- append(class(DOLS_k), "workflow")

  DOLS_k

}

Example Case

Say we have the cointegrating relationship: MB ~ RTPS from lmtest::valueofstocks. (The interpretation is not important here).

dols <- buildDOLS(coint_formula = MB ~ RTPS, data = valueofstocks, fixedk = NULL, robusterrors = T, selection = BIC)

dols 

summary(dols) # the Non-HAC, biased standard errors but model fit results

dols$selection # shows the model selection process

dols$k # the result of that process

dols$call # a stand-alone call function to replicate results with dynlm

dols$HAC # the HAC estimated standard errors/significance values 

Output

The output is not ideal as it is hidden within the model. I have written a few methods but would like some feedback. Does it clearly show what is interesting to the user? Is it too long? etc.

Print Method Shows model selection and a cleaner eval

print.workflow <- function(x) {

  if(!is.null(x$selection)){

    cat("n    Selecting k for Modeln")

    cat("______________________________n")

    print(x$selection)

  }

  cat("n    Model with k =",x$k,"n")

  cat("______________________________")

  print(eval(x$call))

}

dols #** Is it too long? Does it give the important information?

Summary Method Shows HAC estimated errors automatically

summary.workflow <- function(x, ...){

  if("HAC" %in% names(x)){ # if the function has a robusterrors arg

    addHAC <- NextMethod(x) # alter the summary coeffs

    addHAC$coefficients <- x$HAC

    cat("*Summary table depicts HAC estimated errors found by:n")

    cat(paste0("lmtest::coeftest(model, vcov = sandwich::NeweyWest(model, lag = ",x$k,"))n"))

    addHAC

  } else {

    warning("tSummary table does not depict HAC estimated errorsntPlease indicate the buildDOLS robusterrors argument to be TRUE")

    print(NextMethod(x))

  }

} 

summary(dols)

Please let me know what you think!

Nice package, looks pretty solid.

The choice of k, half the cube root, is not obvious. Please cite http://www.ssc.wisc.edu/~kwest/publications/1990/ Automatic Lag Selection eqn. 2.1 (though I still found gamma hat & parameter a bit mysterious).

The "maximum number of lags/leads" comment is helpful. Could we please promote it to a @param appearing near the function signature, so Roxygen2 / devtools::document() will find it?