View source: R/family.others.R

AR1EIM | R Documentation |

Computation of the exact Expected Information Matrix of
the Autoregressive process of order-`1`

(AR(`1`

))
with Gaussian white noise and stationary
random components.

```
AR1EIM(x = NULL, var.arg = NULL, p.drift = NULL,
WNsd = NULL, ARcoeff1 = NULL, eps.porat = 1e-2)
```

`x` |
A vector of quantiles. The gaussian time series for which the EIMs are computed. If multiple time series are being analyzed, then |

`var.arg` |
Logical. Same as with |

`p.drift` |
A numeric vector with the |

`WNsd, ARcoeff1` |
Matrices.
The standard deviation of the white noise, and the
correlation (coefficient) of the AR( That is, the dimension for each matrix is |

`eps.porat` |
A very small positive number to test whether the standar deviation
( See below for further details. |

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`

.

An array of dimension `N \times NOS \times 6`

,
as above.

This array stores the EIMs calculated from the joint density as a function of

`\boldsymbol{\theta} = (\mu^*, \sigma^2, \rho). `

Nevertheless, note that, under the VGAM/VGLM approach, the EIMs
must be correspondingly calculated in terms of the linear
predictors, `\boldsymbol{\eta}`

.

For large enough `n`

, the EIMs,
`J_n(\boldsymbol{\theta})`

,
become approximately linear in `n`

. That is, for some
`n_0`

,

```
J_n(\boldsymbol{\theta}) \equiv
J_{n_0}(\boldsymbol{\theta}) + (n - n_0)
\bar{J}(\boldsymbol{\theta}),~~~~~~(**)
```

where ` \bar{J}(\boldsymbol{\theta}) `

is
a constant matrix.

This relationsihip is internally considered if a proper value
of `n_0`

is determined. Different ways can be adopted to
find `n_0`

. In `AR1EIM`

, this is done by checking
the difference between the internally estimated variances and the
entered ones at `WNsd`

.
If this difference is less than
`eps.porat`

at some iteration, say at iteration `n_0`

,
then `AR1EIM`

takes
` \bar{J}(\boldsymbol{\theta})`

as the last computed increment of
`J_n(\boldsymbol{\theta})`

, and extraplotates
`J_k(\boldsymbol{\theta})`

, for all
`k \geq n_0 `

using `(*)`

.
Else, the algorithm will complete the iterations for
`1 \leq n \leq N`

.

Finally, note that the rate of convergence reasonably decreases if
the asymptotic relationship `(*)`

is used to compute
`J_k(\boldsymbol{\theta})`

,
`k \geq n_0 `

. Normally, the number
of operations involved on this algorithm is proportional to
`N^2`

.

See Porat and Friedlander (1986) for full details on the asymptotic behaviour of the algorithm.

Arguments `WNsd`

, and `ARcoeff1`

are matrices of dimension
`N \times NOS`

. Else, these arguments are accordingly
recycled.

For simplicity, one can assume that the time series analyzed has
a 0-mean. Consequently, where the family function
`AR1`

calls `AR1EIM`

to compute
the EIMs, the argument `p.drift`

is internally set
to zero-vector, whereas `x`

is *centered* by
subtracting its mean value.

V. Miranda and T. W. Yee.

Porat, B. and Friedlander, B. (1986).
Computation of the Exact Information Matrix of Gaussian Time Series
with Stationary Random Components.
*IEEE Transactions on Acoustics, Speech, and Signal Processing*,
**54(1)**, 118–130.

`AR1`

.

```
set.seed(1)
nn <- 500
ARcoeff1 <- c(0.3, 0.25) # Will be recycled.
WNsd <- c(exp(1), exp(1.5)) # Will be recycled.
p.drift <- c(0, 0) # Zero-mean gaussian time series.
### Generate two (zero-mean) AR(1) processes ###
ts1 <- p.drift[1]/(1 - ARcoeff1[1]) +
arima.sim(model = list(ar = ARcoeff1[1]), n = nn,
sd = WNsd[1])
ts2 <- p.drift[2]/(1 - ARcoeff1[2]) +
arima.sim(model = list(ar = ARcoeff1[2]), n = nn,
sd = WNsd[2])
ARdata <- matrix(cbind(ts1, ts2), ncol = 2)
### Compute the exact EIMs: TWO responses. ###
ExactEIM <- AR1EIM(x = ARdata, var.arg = FALSE, p.drift = p.drift,
WNsd = WNsd, ARcoeff1 = ARcoeff1)
### For response 1:
head(ExactEIM[, 1 ,]) # NOTICE THAT THIS IS A (nn x 6) MATRIX!
### For response 2:
head(ExactEIM[, 2 ,]) # NOTICE THAT THIS IS A (nn x 6) MATRIX!
```

VGAM documentation built on Sept. 19, 2023, 9:06 a.m.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.