Implementing VAR predictions manually I have trained a VAR model via statsmodels library in Python, and I am trying to implement the one-step ahead prediction manually without using the package.
Based on statsmodels documentation here, they define a VARMAX process as:
$$y_t = A(t) + A_1 y_{t-1} + \dots + A_p y_{t-p} + B x_t + \epsilon_t +
M_1 \epsilon_{t-1} + \dots M_q \epsilon_{t-q}$$
In my case I do not have exogenous variables nor the MA terms so the model should simplify to a VAR model:
$$y_t = A(t) + A_1 y_{t-1} + \dots + A_p y_{t-p} + \epsilon_t$$
My sketch of Python implementation for prediction is as follows for an example VAR(3) model:
values = []
for i in range(5000):
    y_t = intercept + np.matmul(A_1, lag_1) + np.matmul(A_2, lag_2) + np.matmul(A_3, lag_3) + np.random.multivariate_normal([0, 0], covariance_matrix).reshape(-1, 1)
    values.append(y_t)

where intercept and $y_t$ are $k \times 1$ vectors, $A_1, A_2, A_3$ are parameter matrices, and $lag_1, lag_2, lag_3$ are vectors of $i$-th previous values. I run the equation 5000 times since I sample from the multivariate normal distribution and take the mean of $y_t$, for one-step ahead prediction. I take all the parameters from the trained model.
My estimates here are close to the one-step ahead prediction statsmodels gives but are not quite the same, so I am not doing something right. Advice is appreciated.

VAR predictions are computed as a conditional expectation, e.g.
$$E[y_{t+h}|y_t, y_{t-1}, \dots] = A(t) + A_1 y_t + \dots + A_p y_{t-p+1}$$
The error term $\epsilon_t$ does not enter into the prediction because it is defined such that $E[\epsilon_{t+h}|y_t, y_{t-1}, \dots] = 0$.
Even with 5000 samples, your sample mean $\frac{1}{n} \sum_{i=1}^{5000} e_t$ will generally not be exactly equal to zero, and so your prediction will differ slightly from the one produced by statsmodels.
