Is there a function like auto.arima which gives the best VAR model according to a metric like AIC? The function auto.arima from the forecast package, automatically fits the best ARIMA model given a time series, by evaluating on a metric like AIC, AICc or BIC.
I am looking for a similar function which does the same, but fits a VAR (Vector Auto Regression) model to a bunch of time series variables. Is there any such function in any package?

A close one is VARselect from the vars package. It returns optimal lag orders according four information criteria (AIC, BIC a.k.a. SC, FPE and HQ). Then you can manually specify the selected lag order and estimate the model with VAR as follows:
aic <- VARselect(Canada, lag.max=5, type="const")$selection["AIC(n)"]
model <- VAR(Canada, p=aic, type="const")
print(model)

