Upon the multivariate representation, $$\Phi_0 y_t = \Psi + \Phi_1 Y_{T-1} + ... + \Phi_P y_{T-P} + \epsilon_T ,$$
where the $\Phi_i, i=1,2,...,P$ are $s \times s$ matrices containing the $\phi_{is}
parameters.$, the one-step-ahead forecasts for the year $T+1$ is straightforward,
$$y_t = \Phi_0^{-1} \Psi + \Phi_0^{-1} \Phi_1 Y_{T-1} + ... + \Phi_0^{-1} \Phi_P y_{T-P} +
\Phi_0^{-1} + \epsilon_T .$$
Multi-step-ahead forecasts are obtained recursively.
The prediction errors variances for the one-step-ahead forecast are the diagonal elements of
$$\sigma^2 \Phi_0^{-1} (\Phi_0^{-1})^{'},$$
whereas for $h=2,3,...$ years ahead forecasts it becomes
$$\sigma^2 \Phi_0^{-1} (\Phi_0^{-1})^{'} + (h-1) (\Gamma \Phi_0^{-1}) (\Gamma \Phi_0^{-1})^{'},$$
where $\Gamma = \Phi_0^{-1} \Phi_1$.
This version considers PIAR models up to order 2 for quarterly observed data. By default, seasonal
intercepts are included in the model as deterministic components.
The number of observations to forecast, hpred must be a multiple of 4.