This function implements the algorithm of Porat and Friedlander
(1986) to recursively compute the exact expected
information matrix (EIM) of Gaussian time series with stationary
random components.
By default, when the VGLM/VGAM family function
AR1 is used to fit an AR(\(1\)) model
via vglm, Fisher scoring is executed using
the approximate EIM for the AR process. However, this model
can also be fitted using the exact EIMs computed by
AR1EIM.
Given \(N\) consecutive data points,
\( {y_{0}, y_{1}, \ldots, y_{N - 1} } \) with probability density \(f(\boldsymbol{y})\),
the Porat and Friedlander algorithm
calculates the EIMs
\( [J_{n-1}(\boldsymbol{\theta})] \),
for all \(1 \leq n \leq N\). This is done based on the
Levinson-Durbin algorithm for computing the orthogonal polynomials of
a Toeplitz matrix.
In particular, for the AR(\(1\)) model, the vector of parameters
to be estimated under the VGAM/VGLM approach is
$$ \boldsymbol{\eta} = (\mu^{*}, \log(\sigma^2), rhobit(\rho)),$$
where \(\sigma^2\) is the variance of the white noise
and \(mu^{*}\) is the drift parameter
(See AR1 for further details on this).
Consequently, for each observation \(n = 1, \ldots, N\), the EIM,
\(J_{n}(\boldsymbol{\theta})\), has dimension
\(3 \times 3\), where the diagonal elements are:
$$ J_{[n, 1, 1]} =
E[ -\partial^2 \log f(\boldsymbol{y}) / \partial ( \mu^{*} )^2 ], $$
$$ J_{[n, 2, 2]} =
E[ -\partial^2 \log f(\boldsymbol{y}) / \partial (\sigma^2)^2 ], $$
and
$$ J_{[n, 3, 3]} =
E[ -\partial^2 \log f(\boldsymbol{y}) / \partial ( \rho )^2 ]. $$
As for the off-diagonal elements, one has the usual entries, i.e.,
$$ J_{[n, 1, 2]} = J_{[n, 2, 1]} =
E[ -\partial^2 \log f(\boldsymbol{y}) / \partial \sigma^2
\partial \rho], $$
etc.
If var.arg = FALSE, then \(\sigma\) instead of \(\sigma^2\)
is estimated. Therefore, \(J_{[n, 2, 2]}\),
\(J_{[n, 1, 2]}\), etc., are correspondingly replaced.
Once these expected values are internally computed, they are returned
in an array of dimension \(N \times 1 \times 6\),
of the form
$$J[, 1, ] = [ J_{[ , 1, 1]}, J_{[ , 2, 2]}, J_{[ , 3, 3]},
J_{[ , 1, 2]}, J_{[, 2, 3]}, J_{[ , 1, 3]} ]. $$
AR1EIM handles multiple time series, say \(NOS\).
If this happens, then it accordingly returns an array of
dimension \(N \times NOS \times 6 \). Here,
\(J[, k, ]\), for \(k = 1, \ldots, NOS\), is a matrix
of dimension \(N \times 6\), which
stores the EIMs for the \(k^{th}\)th response, as
above, i.e.,
$$J[, k, ] = [ J_{[ , 1, 1]}, J_{[ , 2, 2]},
J_{[ , 3, 3]}, \ldots ], $$
the bandwith form, as per required by
AR1.