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
# **determine a maximum lead or lag to start with**.
# 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
k <- floor(nrow(M1Germany)^(1/3))/2 # max lags (or leads) to start with.
# 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
# [1] 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!
