longmemo (version 1.1-3)

WhittleEst: Whittle Estimator for Fractional Gaussian Noise / Fractional ARIMA

Description

Computes Whittle's approximate MLE for fractional Gaussian noise or fractional ARIMA (=: fARIMA) models, according to Beran's prescript.

Relies on minmizing Qeta() (\(= \tilde{Q}(\eta)\), which itself is based on the “true” spectrum of the corresponding process; for the spectrum, either specFGN or specARIMA is used.

Usage

WhittleEst(x,
           periodogr.x = per(if(scale) x/sd(x) else x)[2:((n+1) %/% 2)],
           n = length(x), scale = FALSE,
           model = c("fGn", "fARIMA"),
           p, q,
           start = list(H= 0.5, AR= numeric(), MA=numeric()),
           verbose = getOption("verbose"))
# S3 method for WhittleEst
print(x, digits = getOption("digits"), ...)

Value

An object of class WhittleEst which is basically a list with components

call

the function call.

model

= input

n

time series length length(x).

p,q

for "fARIMA": order of AR and MA parts, respectively.

coefficients

numeric 4-column matrix of coefficients with estimate of the full parameter vector \(\eta\), its standard error estimates, z- and P-values. This includes the Hurst parameter \(H\).

theta1

the scale parameter \(\hat{\theta_1}\), see Qeta.

vcov

the variance-covariance matrix for \(\eta\).

periodogr.x

= input (with default).

spec

the spectral estimate \(\hat{f}(\omega_j)\).

There is a print method, and coef,

confint or vcov methods work as well for objects of class "WhittleEst".

Arguments

x

numeric vector representing a time series. Maybe omitted if periodogr.x and n are specified instead.

periodogr.x

the (raw) periodogram of x; the default, as by Beran, uses per, but tapering etc may be an alternative, see also spec.pgram.

n

length of the time series, length(x).

scale

logical indicating if x should be standardized to (sd) scale 1; originally, scale = TRUE used to be built-in; for compatibility with other methods, notably plotting spectra, scale = FALSE seems a more natural default.

model

numeric vector representing a time series.

p,q

optional integers specifying the AR and MA orders of the fARIMA model, i.e., only applicable when model is "fARIMA".

start

list of starting values; currently necessary for model = "fARIMA" and with a reasonable default for model = "fGn".

verbose

logical indicating if iteration output should be printed.

digits,...

optional arguments for print method, see print.default.

Author

Martin Maechler, based on Beran's “main program” in Beran(1994).

References

Beran, Jan (1994). Statistics for Long-Memory Processes; Chapman & Hall. (Section 6.1, p.116--119; 12.1.3, p.223 ff)

See Also

Qeta is the function minimized by these Whittle estimators.

FEXPest for an alternative model with Hurst parameter, also estimated by a “Whittle” approximate MLE, i.e., a Whittle's estimator in the more general sense.

The plot method, plot.WhittleEst.

Examples

data(NileMin)
(f.Gn.N  <- WhittleEst(NileMin))                             # H = 0.837
(f.A00.N <- WhittleEst(NileMin, model = "fARIMA", p=0, q=0)) # H = 0.899
confint(f.Gn.N)
confint(f.A00.N)

data(videoVBR)
(f.GN    <- WhittleEst(videoVBR))

## similar {but faster !}:
(f.am00  <- WhittleEst(videoVBR, model = "fARIMA", p=0, q=0))
rbind(f.am00$coef,
      f.GN  $coef)# really similar

f.am11  <- WhittleEst(videoVBR, model = "fARIMA",
                      start= list(H= .5, AR = .5, MA= .5))
f.am11
vcov(f.am11)

op <- if(require("sfsmisc"))
  mult.fig(3, main = "Whittle Estimators for videoVBR data")$old.par  else
  par(mar = c(3,1), mgp = c(1.5, 0.6, 0), mar = c(4,4,2,1)+.1)
plot(f.GN)
plot(f.am00)
plot(f.am11)

et <- as.list(coef(f.am11))
et$AR <- c(et$AR, 0, 0) # two more AR coefficients ..
f.am31  <- WhittleEst(videoVBR, model = "fARIMA", start = et)
## ... warning non nonconvergence,  but "kind of okay":
lines(f.am31, col = "red3") ## drawing on top of ARMA(1,1) above - *small* diff

f.am31 # not all three are "significant"
round(cov2cor(vcov(f.am31)), 3) # and they are highly correlated

et <- as.list(coef(f.am31))
et$AR <- unname(unlist(et[c("AR1", "AR2")]))
f.am21  <- WhittleEst(videoVBR, model = "fARIMA",
                      start = c(et[c("H","AR", "MA")]))
f.am21
lines(f.am21, col = adjustcolor("gold", .4), lwd=4)

par(op)## (reset graphic layout)

##--> ?plot.WhittleEst  for an example using  'periodogr.x'