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.