Estimating VARX(p,q) models in Mathematica

I would like to estimate the parameters of a vector autogressive process VARX(p,q) with n endogenous variables y and m exogenous variables x: $$y_t = v + B_1 y_{t-1} + ... + B_p y_{t-p} + \theta_0 x_t + ... + \theta_q x_{t-q} + e_t$$ where $$e_t \sim N(0,\Sigma)$$

The Mathematica functions TimeSeriesModelFit[data,mspec] and EstimatedProcess[data,proc] only accept ARMA, ARIMA and GARCH processes without exogenous variables.

What would be the best approach to estimate the coefficients $B_i$ and $\theta_i$? We could transform the system as $Y=B \cdot Z + U$ where $B = [v\space B_1\space B_2 ... B_p\space \theta_0\space \theta_1 ... \theta_q]$ and solve it by OLS. Is there a way to use the existing functions TimeSeriesModelFit and EstimatedProcess for processes with exogenous variables?