# Breusch-Godfrey LM Test returning different results in vars and lmtest

I like to test a time series for serial correlation in order to determine the optimum lag length of a VAR. I created following reproducible example:

library(vars)
library(lmtest)

# Data

# Optimum lag length (lag.max = 4 because we have quarterly data)
VARselect(Canada[,2:3], lag.max = 4, type = "const")\$selection
# 1 or 3 ?

# Test serial correlation in the residuals with the Breusch-Godfrey LM test


serial.test from vars

serial.test(VAR(Canada[,2:3], p = 1, type = "const"), type="BG")

Breusch-Godfrey LM test
Chi-squared = 33.164, df = 20, p-value = 0.03237

serial.test(VAR(Canada[,2:3], p = 2, type = "const"), type="BG")

Breusch-Godfrey LM test
Chi-squared = 22.456, df = 20, p-value = 0.3163

serial.test(VAR(Canada[,2:3], p = 3, type = "const"), type="BG") # 3?

Breusch-Godfrey LM test
Chi-squared = 24.313, df = 20, p-value = 0.229

serial.test(VAR(Canada[,2:3], p = 4, type = "const"), type="BG")

Breusch-Godfrey LM test
Chi-squared = 30.27, df = 20, p-value = 0.06559


bgtest from lmtest

bgtest(residuals(VAR(Canada[,2:3], p = 1, type = "const"))[,1] ~ residuals(VAR(Canada[,2:3], p = 1, type = "const"))[,2])

Breusch-Godfrey test for serial correlation of order up to 1
LM test = 5.7701, df = 1, p-value = 0.0163

bgtest(residuals(VAR(Canada[,2:3], p = 2, type = "const"))[,1] ~ residuals(VAR(Canada[,2:3], p = 2, type = "const"))[,2])

Breusch-Godfrey test for serial correlation of order up to 1
LM test = 0.1004, df = 1, p-value = 0.7514

bgtest(residuals(VAR(Canada[,2:3], p = 3, type = "const"))[,1] ~ residuals(VAR(Canada[,2:3], p = 3, type = "const"))[,2]) # 3?

Breusch-Godfrey test for serial correlation of order up to 1
LM test = 0.022743, df = 1, p-value = 0.8801

bgtest(residuals(VAR(Canada[,2:3], p = 4, type = "const"))[,1] ~ residuals(VAR(Canada[,2:3], p = 4, type = "const"))[,2])

Breusch-Godfrey test for serial correlation of order up to 1
LM test = 0.073025, df = 1, p-value = 0.787


Questions:

1. Why does serial.test create different results than bgtest?
2. Would you agree that lag = 3 is the best choice?

You are submitting univariate series to bgtest and multivariate series to serial.test. Since the data being tested differs, why should the test results be the same? They should not, and they are not.
You can also see from the output that the degrees of freedom differ between bgtest and serial.test. That immediately suggests that the test specifications differ.
Which lag is the best choice might not be clear from the test results. You can always increase the lag of the VAR model so much that you severely overfit and the residuals become roughly uncorrelated. But overfitting is not desirable. You want not only good performance on the training data, but also out of sample. You could look at information criteria to pick the "best" lag, as facilitated by the function VARselect in the "vars" package. This is not a magic bullet, but at least it offers a balance between overfitting and underfitting.