# [R]: DOLS - number of leads and lags

I have a model

model_1=dynlm(y_log~x1_log+x2)


where y_log, x1_log are I(1) and x2 is something between I(0) and I(1) (I rejected unit root with ADF test but KPSS test rejected stationarity).

I would like to estimate this regression and get coefficients I could trust. I am looking at this regression as on cointegration (from books I know that to test for cointegration relationship using Engle-Granger procedure all variables should be of the same order but I was advised that this rule is not so strict when I have at least 2 variables that are I(1)). Residuals from this regression are stationary that implies cointegration.

I was advised that in my case it would be appropriate to run a dynamic OLS, i.e. to add some leads and lags of the differences of variables that are I(1) - based on Stock and Watson (1993). However, I am totally lost on how should I choose the number of leads and lags. I read that this could be done based on AIC but I don't entirely understand how. Should I run various models with different number or leads and lags and then to compare them or is there a command in R that would give me the number straight away?

When I am adding leads and lags, the AIC value is rising and the equation without any leads and lags has the smallest AIC. However, if I am adding just the lagged values, the AIC is falling and I can't find the lag at which it stops falling.

Could someone, please, advise me how I should proceed and what functions in R I might use? Thanks.

Here is an example of what I typically do, in code form easy running:

# Loading Packages and Data
library(dynlm)
data("M1Germany")
# convert to time series object
M1Germany <- window(as.ts(M1Germany, start = c(1960,1), frequency = 4), end = c(1995,4))
# A basic OLS equation with nonstationary data assuming you have found they cointigrate
# - here they do not cointigrate but I wanted to use an easily loaded dataset~
dynlm(loggnp ~ logprice + interest, data = M1Germany)
# We can do better.

#There are a few steps to think about:

# **1.**
# In order to test a few different lags/leads, you must
# Every statistician seems to have a different ad hoc method of doing this--
# I learned from someone to just use the
# cubed root of the number of observations as a starting point, divided by 2

# 2 to start (2 lags and 2 leads)

# Typically when I use a model that only has lagged values...
# I use AIC (like for an adf test), and only floor(nrow(M1Germany)^(1/3)) , no "/2"
# **but when I incorporate lead values (DOLS model)...
# I use BIC** - also referred to as the Schwarz Bayesian Criterion (SBC).
# You will want to choose the model that minimizes the BIC value.

bic <- function (model) { # what I use in my own package
n <- df.residual(model) + length(variable.names(model))
log(sum(resid(model)^2)/n) + length(variable.names(model)) * log(n)/n
}

# **2.** I
#mportantly, you will want to **fix the dates of the models
# you are using to the same number of observations**.
# Here I use the dynlm package and start().

DOLS <- dynlm(loggnp ~ logprice + interest +
L(diff(logprice),-k:k) +
L(diff(interest),-k:k),
data = M1Germany) # a DOLS model that builds lags and leads
start(DOLS) # the date to start, 1961 March
#  1961    3

#Then, in order to test a bunch of different lags on your model
# you can go ahead and use something along the lines of:
bics <- sapply(1:k, function(k) bic(dynlm(loggnp ~ logprice + interest +
L(diff(logprice),-k:k) +
L(diff(interest),-k:k),
data = M1Germany, start = start(DOLS))))
# Where you cycle through using 1:k as our model's -k:k s
# and fix the start date with start = start(DOLS)
k <- which.min(bics)  # it's 1 lag/lead

summary(dynlm(loggnp ~ logprice + interest +
L(diff(logprice),-k:k) +
L(diff(interest),-k:k),
data = M1Germany))
# Time series regression with "ts" data:
#   Start = 1960(3), End = 1995(3)
#
# Call:
#   dynlm(formula = loggnp ~ logprice + interest + L(diff(logprice),
#                                                    -k:k) + L(diff(interest), -k:k), data = M1Germany)
#
# Residuals:
#   Min        1Q    Median        3Q       Max
# -0.160804 -0.040456  0.009372  0.035920  0.091552
#
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
#   (Intercept)                6.53504    0.05111 127.871  < 2e-16 ***
#   logprice                   0.55475    0.01169  47.435  < 2e-16 ***
#   interest                   0.49022    0.38543   1.272 0.205654
#   L(diff(logprice), -k:k)-1  1.08651    0.27305   3.979 0.000114 ***
#   L(diff(logprice), -k:k)0   1.16442    0.30016   3.879 0.000165 ***
#   L(diff(logprice), -k:k)1   0.41560    0.26366   1.576 0.117351
#   L(diff(interest), -k:k)-1  2.10670    0.97471   2.161 0.032471 *
#   L(diff(interest), -k:k)0   1.69783    1.00642   1.687 0.093966 .
#   L(diff(interest), -k:k)1   1.68935    0.99369   1.700 0.091473 .
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.05424 on 132 degrees of freedom
# Multiple R-squared:  0.9482,  Adjusted R-squared:  0.945
# F-statistic: 301.8 on 8 and 132 DF,  p-value: < 2.2e-16


Remember however that the standard error (and thus the t test and p value too) are unreliable and need to be estimated using robust errors techniques - refer to the coeftest in the sandwich package.

The whole process can be a little tedious. I myself have been writing my own package for these kinds of exploratory time series econometrics methods. If you are interested it can be found on GitHub constantly being edited as I learn more about econometrics.

install.packages("remotes")
remotes::install_github("efriedland/friedland")
library(friedland)
data("MacKinnon")
buildDOLS(loggnp ~ logprice + interest, data = M1Germany)


Any questions feel free to ask!

• Thank you very much for you answer! It was very helpful. Cheers – Marika May 9 '17 at 7:59
• Happy to help. If you found this answer satisfactory feel free to mark the check symbol as answered~ – Red May 9 '17 at 12:29