#You can just run the code in your console
library(quantmod)
library(vars)


#Download the two stock symbols in question VNQ and SPY
tickers = c("SPY", "VNQ")
getSymbols(tickers)

#Convert data to a stationairy process
spy_station = diff(SPY$SPY.Close)
vnq_station = diff(VNQ$VNQ.Close)

#merge the data together
data = merge(spy_station, vnq_station)
data = na.omit(data)

#Now we see which lag is optimal
VARselect(data, type = "both")
    1         2         3         4         5         6         7         8
    9        10
    AIC(n) 0.4611156 0.4602640 0.4605474 0.4611319 0.4614260 0.4630208
    0.4637955 0.4607941 0.4623594 0.4616697
    #It looks like p = 2 is the best.

 #Create our VAR model
 vmodel = VAR(data, p = 2, type = "both")
 summary(vmodel)