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:


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
    # 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")

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)
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 


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) {
    cat("\n    Selecting k for Model\n")
  cat("\n    Model with k =",x$k,"\n")
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"))
  } else {
    warning("\tSummary table does not depict HAC estimated errors\n\tPlease indicate the buildDOLS robusterrors argument to be TRUE")

Please let me know what you think!


1 Answer 1


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?

  • \$\begingroup\$ Hi J H thanks for looking. I I wonder if you are talking about the package I put this in or just the above code? Also, thank you for tracking down a citation (I myself have not read this paper). Assuming it provides an optimal number of maximum lags to start selecting from, because this is a lag AND lead model, the cube root would tell me 4 lags, when the model is really looking for 2 lags 2 leads (k). Does this make sense? And where would you propose promoting*? I currently have not adapted my package to be built by Roxygen2, only manually changing items on Github (efriedland/friedland) \$\endgroup\$ Sep 10, 2017 at 21:23
  • \$\begingroup\$ I was just looking at the posted code. The paper is not completely transparent, but it does make it clear that the desired exponent will change with kernel type, and 1/3 (or 2/9) is specifically for Bartlett. It sounds like you got the 0.5 coefficient from elsewhere. I suspect 2 + 2 leads / lags is a lot like 4 lags, if there's enough data that we're not worried about the edge cases. By promote I meant list at top where Help can find it, similar to javadoc or python docstrings. Basically I was calling it a useful comment, worthy of being more easily accessed. \$\endgroup\$
    – J_H
    Sep 10, 2017 at 21:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.