Multivariate Time Series Analysis Approach

I am facing an issue regarding the overall approach required for multivariate time series forecasting. I am a novice in R and statistics. Suppose I have 3 time series, X, Y and Z, where Z depends on either X alone or Y alone or difference of X and Y. My objective is to forecast Z using X or Y or X-Y. Now from what I understood from various blogs and online materials, first perform a granger causality check to find whether Z can be better forecasted with past values of Z,X and Y. If so, perform VAR, else check the cointegration between the time series using Johansen test and if there is a cointegration ( which suggests causality ), perform VECM and then predict it. If the cointegration test suggest no causation, then perform normal Arima with external regressors or Arima with regression errors (as mentioned in Otext forecasting principles) and check the forecast accuracy.

Now below are my queries.

1) Is the approach I have mentioned above correct, if not, kindly outline the correct approach

2) Is granger causality test a prerequisite for VAR, because from what I read, VAR is used whenever both or multiple time series influence each other. If gc test is significant for Y~X and not vice versa, can I proceed with VAR? Also, it sounds like gc test is more like a prerequisite for transfer functions.

3) In granger causality, if X is stationary whereas Y is not , should I apply differencing only to Y or use the same order differencing to both X and Y prior to the test.

4) For finding transfer function in r, I have used arfima in the below code,

data(SeriesJ)
YJ <- SeriesJ$YJ XJ <- SeriesJ$XJ
prewhiten(XJ,YJ)
fitTF1 <- arfima(YJ, order= c(2, 0, 0), xreg = XJ, reglist = list(regpar = c(2, 2, 3)), lmodel = "n")

where regpar depends on the cross correlation plot between xj and yj.However, this is giving me an error. Kindly suggest an alternative command for doing transfer functions. I read that R has built in ARIMAX called arima, so will this command suffice

fitTF1 <- arima (Y, order=c(1, 2, 3), seasonal=list (order=c(1, 0, 0),     freq=12), xreg=X)
pred <- predict (fitTF1,9,newxreg=X.predict)

5) Suppose, the above queries are resolved and I found Arima with regression errors to be a better multivariate model than VAR and VECM and I have decided to go ahead with it. So,I initially split the data ( X, Y and Z) intro training and testing set,build the model using training set, validate using test set. Now for actual forecasting, I am considering entire, X with Y and Z as regressors and then forecast. I have actual values for Y in the forecast period(Ex post forecast) and I don’t have future values for Z, hence I need to do univariate forecasting for Z and use the forecasted values as regressors. Is the approach I outlined correct, moreover, say the model tested well in the test set with low RMSE and MAPE value, but the forecast of Z is coming very poorly ( Z has intermittent data or no seasonality) and accuracy metrics in test set of Z are giving poor values, in such a case, should we continue using forecast of Z, or stick with Y alone as a regressor.

Hope you understood my queries. Kindly let me know if you need more clarity.