I have some trouble understanding the forecasting/inference process of ARMA models.
From Hamilton (which I am reading now), we can obtain forecasts at $Y$ from any linear process with r.v. values $X$ using the expectations as follows:
$\mathbf{a} = E(XX^T)^{-1}E(XY^T)$
This assumes that the forecast is the optimal linear forecast using projections. The matrix $E(XX^T)$ is the matrix of autocovariances (since ARMA is stationary), and $E(XY^T)$ is also a matrix of autocovariances.
1) We can establish all autocovariances from the sample (time series) and forecast the future values using the above formula by simply substituting estimates for second moments, inverting matrices, and we are done. Why do we need to find the ARMA coefficients to forecast which is done in all packages like stat::arima?
2) I still cannot figure the whole process - stat::arima uses the Kalman filter to compute the fundamental innovations, and somehow compute the likelihood (cannot decipher Hamilton chapter 13 yet since have not yet read the prerequisite chapters). From stat::arima:
The exact likelihood is computed via a state-space representation of the ARIMA process, and the innovations and their variance found by a Kalman filter.
So, these innovations + observations should be enough to forecast the future values. Why do we even need likelihood at all? From Gardner et al. https://www.jstor.org/stable/2346910:
The prediction and updating are carried out by means of a set of recursive equations known as the "Kalman filter" The parameter $\sigma^2$ does not appear in the recursion
3) When Kalman filter finds the innovations - it finds also the variance of those and computes the maximum likelihood value. BUT HOW ARE PARAMETERS of ARMA for maximum likelihood found? Kalman can compute the maximum likelihood values, but how are the maximum likelihood parameters $\theta$ and $\phi$ selected? As far as I know, to find these parameters, some gradient descent has to be used, or Newton method. From Gardner again:
The log-likelihood function may then be maximized with respect to ($\theta, \phi$) by minimizing:
$\mathcal{L}(\theta, \phi) = n \log S(\theta, \phi) + \sum_{i = 1}^n \log f_i, \quad \quad f_i \propto \hat{MSE}$
So, the procedure described in the paper outputs the values of errors. likelihood, etc. But where does it get the optimal $\theta, \phi$?