# What is the proper procedure for forecasting a VAR estimated in differences?

I am comparing the forecasting accuracy of a VAR estimated in levels and first differences. I am working in R, but my question is language agnostic (although if you have an answer with R code, that will be appreciated.)

Specifically, when estimating a VAR in levels, the forecast that you'll get from your statistical package will be in the levels of that variable. If you care about the future levels of a time series, that is great, you're done.

If you forecast in first-differences, the statistical package will likely provide you forecasts of future first-differences in the time series, not the levels. So my question is if I estimate a VAR in first-differences, is there any problem with just applying the forecasted first-differences iteratively with an initial value equal to the last observed data point in the time-series? Is this equivalent to taking the first-differenced VAR model (i.e. the coefficients) and applying it to the level data to obtain forecasted levels?

For R users, I'm asking if I can estimate a VAR in first differences, call it var_fd, and pass it to the forecast function, and then back out the forecasted levels using diffinv and my last observed data point.

No problem there. For an $h$-step-ahead forecast you may simply take the $h$ predicted differences for time periods $t+1,\dots,t+h$ and add them to the last observed value $y_t$ to obtain the forecast $\hat y_{t+h|t}$ for $y_{t+h}$ as of time $t$: $$\hat y_{t+h|t} = y_t + \widehat{\Delta y}_{t+1|t} + \dotsc + \widehat{\Delta y}_{t+h|t}$$
No. Here is a counterexample: consider an estimated model in first differences, $$\begin{pmatrix} \Delta y_{1,t} \\ \Delta y_{2,t} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} \Delta y_{1,t-1} \\ \Delta y_{2,t-1} \\ \end{pmatrix} + \begin{pmatrix} v_{1,t} \\ v_{2,t} \\ \end{pmatrix}.$$ It will produce forecasts $\widehat{\Delta y}_{i,t+h|t}=h\Delta y_{i,t}$ and thus $\hat y_{i,t+h|t}=y_{i,t}+h\Delta y_{i,t}$ for $i=1,2$.
However, when the estimated coefficients are applied on data in levels, $$\begin{pmatrix} y_{1,t} \\ y_{2,t} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix} \begin{pmatrix} y_{1,t-1} \\ y_{2,t-1} \\ \end{pmatrix} + \begin{pmatrix} u_{1,t} \\ u_{2,t} \\ \end{pmatrix},$$ the forecasts are $\hat y_{i,t+h|t}=y_{i,t}$ for $i=1,2$. The two are not the same.