The development of advanced recording and measurement devices in scientific fields are producing high-dimensional time series data. Vector autoregressive (VAR) models are well suited for analyzing the high dimensional time series datasets, but accurate inference at scale remains a central challenge. We have recently introduced Union of Intersections (UoI) framework, which in linear models (e.g., UoILasso) exhibits low false-positive and false-negative feature selection along with low bias and variance estimation. Here, we briefly summarize past results, and then extend them in directions required to infer dynamic, causal networks from large quantities of observational data. Through extensive numerical investigation in the context of correlated regressors under varied design matrices, we find that in the sparse regime, UoILasso exhibits favorable feature selection and estimation properties compared to other methods. We next scale the UoILasso algorithm for VAR models to analyze large problem sizes (TBs). These advances enable us to estimate the largest VAR model known (1000 nodes) and apply it to large time-series data from neurophysiology (192 neurons) and finance (470 companies).