Nothing
## Do not edit this file manually.
## It has been automatically generated from *.org sources.
setClass("VirtualAutocovariances", contains = "VIRTUAL")
setClass("VirtualAutocorrelations", contains = "VIRTUAL")
setClass("Autocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("PartialAutocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("Autocorrelations", contains = c("FlexibleLagged", "VirtualAutocovariances"))
setClass("PartialAutocorrelations", contains = c("FlexibleLagged", "VirtualAutocovariances"))
## Commenting out the slots, their intention is valid in the univariate case
setClass("PartialVariances", contains = c("FlexibleLagged", "VirtualAutocovariances")
## , slots = c(data = "numeric")
)
## Combo variants
setClass("ComboAutocovariances", contains = c("FlexibleLagged", "VirtualAutocovariances")
## , slots = c(data = "matrix")
)
setClass("ComboAutocorrelations", contains = c("FlexibleLagged", "VirtualAutocorrelations")
## , slots = c(data = "matrix")
)
setClass("Fitted",
slots = c(n = "numeric", varnames = "character", objectname = "character" )
)
setClass("SampleAutocovariances", contains = c("Autocovariances", "Fitted"))
setClass("SamplePartialAutocovariances", contains = c("PartialAutocovariances", "Fitted"))
setClass("SampleAutocorrelations", contains = c("Autocorrelations", "Fitted"))
setClass("SamplePartialAutocorrelations", contains = c("PartialAutocorrelations","Fitted"))
setClass("SamplePartialVariances", contains = c("PartialVariances", "Fitted"))
## TODO: FittetComboXX variants?
.printFittedPart <- function(x){
for(slot in slotNames("Fitted")){
switch(slot,
varnames = {
if(length(slot(x, slot)) == 0){
cat('Slot ', slot, ': ', " <not set>", "\n", sep = "")
}else{
cat('Slot ', slot, ':', "\n", sep = "")
print(slot(x, slot))
}
},
objectname = {
cat('Slot ', slot, ': ', slot(x, slot), "\n", sep = "")
},
## default
{ cat('Slot ', slot, ':', "\n", sep = "")
print(slot(x, slot))
}
)
}
}
setMethod("show", signature(object = "Autocovariances"),
function (object)
{
.reportClassName(object, "Autocovariances")
# 2019-03-29 (similarly below) was: callNextMethod()
# is it better to drop prefix "Lag_":
# dataWithLagNames(object, "")
print(dataWithLagNames(object))
})
setMethod("show", signature(object = "Autocorrelations"),
function (object)
{
.reportClassName(object, "Autocorrelations")
print(dataWithLagNames(object))
}
)
setMethod("show", signature(object = "PartialAutocovariances"),
function (object)
{
.reportClassName(object, "PartialAutocovariances")
print(dataWithLagNames(object))
}
)
setMethod("show", signature(object = "PartialAutocorrelations"),
function (object)
{
.reportClassName(object, "PartialAutocorrelations")
print(dataWithLagNames(object))
}
)
setMethod("show", signature(object = "PartialVariances"),
function (object)
{
.reportClassName(object, "PartialVariances")
print(dataWithLagNames(object))
}
)
setMethod("show", signature(object = "SampleAutocorrelations"),
function (object)
{
.reportClassName(object, "SampleAutocorrelations")
callNextMethod()
.printFittedPart(object)
}
)
setMethod("show", signature(object = "SampleAutocovariances"),
function (object)
{
.reportClassName(object, "SampleAutocovariances")
callNextMethod()
.printFittedPart(object)
}
)
setMethod("show", signature(object = "SamplePartialAutocovariances"),
function (object)
{
.reportClassName(object, "SamplePartialAutocovariances")
callNextMethod()
.printFittedPart(object)
}
)
setMethod("show", signature(object = "SamplePartialAutocorrelations"),
function (object)
{
.reportClassName(object, "SamplePartialAutocorrelations")
callNextMethod()
.printFittedPart(object)
}
)
setMethod("show", signature(object = "SamplePartialVariances"),
function (object)
{
.reportClassName(object, "SamplePartialVariances")
callNextMethod()
.printFittedPart(object)
}
)
.basecfclass2S4 <- function(x, r0){
stopifnot(inherits(x, "acf")) # 2020-02-29 was: stopifnot(class(x) == "acf"), similarly below
lagged <- Lagged(x) # Lagged knows how to interpret "acf" objects
if(!missing(r0))
lagged[0] <- if(inherits(r0, "acf")) r0$acf[1, , ] else r0
data <- lagged[]
n <- x$n.used
series <- x$series
varnames <- if(is.null(x$snames)) character(0) else x$snames #TODO: NA in the former?
## lag <- x$lag
class <- switch(x$type, covariance = "SampleAutocovariances",
correlation = "SampleAutocorrelations",
partial = "SamplePartialAutocorrelations" )
new(class, data = data, n = n, varnames = varnames, objectname = series)
}
autocovariances <- function(x, maxlag, ...){
if(is.list(x) && all(names(x) %in% c("ar", "ma", "sigma2"))){ # theoretical ARMA autocov
# ARMAacf(x$ar, x$ma, lag.max = maxlag)
# TODO: if using ARMAacf, need to calculate R(0) and multiply by it!
# using the function from package ltsa instead - but need checking;
# does it work for non-invertible models?
## TODO: currently only univariate
phi <- if(is.null(x$ar )) numeric(0) else x$ar
theta <- if(is.null(x$ma )) numeric(0) else x$ma
sigma2 <- if(is.null(x$sigma2)) NA_real_ else x$sigma2
# !!! Note the '-' below!
wrk <- tacvfARMA(phi = phi, theta = - theta, sigma2 = sigma2, maxLag = maxlag)
new("Autocovariances", data = wrk)
}else{ ## numeric, "ts", "mts" etc. ...
# TODO: check what happens if argument 'plot' is in `...' (error?)
wrk <- if(missing(maxlag))
acf(x, plot = FALSE, type = "covariance", ...)
else
acf(x, lag.max = maxlag, plot = FALSE, type = "covariance", ...)
obj <- .basecfclass2S4(wrk)
# TODO: FUNCTIONS IN PCTS PACKAGE WHICH RELY ON RESULT OF TYPE slMatrix()
# IN THE MULTIVARIATE CASE;
# THEY SHOULD CONTERT TO IT THEMSELVES.
#
# res <- acfbase2sl(res)
# slMatrix(res)
obj
}
}
setGeneric("autocovariances")
setGeneric("autocorrelations",
function(x, maxlag, lag_0, ...) standardGeneric("autocorrelations")
)
partialAutocorrelations <- function(x, maxlag, lag_0 = TRUE, ...){
# TODO: needs more work!
# see handwritten notes "parcor"; James Proberts may wish to do it.
#
# in this case, compute autocovariances, solve YW and extract the parcor.
# (for the multivariate case needs modification!)
#
# calculate autocor not autocov, since the latter may need more, e.g. sigma2
acrobj <- if(missing(maxlag))
autocorrelations(x, ...)
else
autocorrelations(x, maxlag = maxlag, ...)
pacr <- modelCoef(acrobj, "PartialAutocorrelations")
if(isTRUE(lag_0)){
if(is(acrobj, "SampleAutocorrelations")){
new("SamplePartialAutocorrelations", data = pacr, as(acrobj, "Fitted"))
}else
new("PartialAutocorrelations", data = pacr)
}else if(identical(lag_0, FALSE)){
# if(dim(res@data) == 3)
# new("FlexibleLagged", data = obj@data[ -1, , ])
# else
# new("FlexibleLagged", data = obj@data[-1])
## TODO: this is very lazy; to avoid checking the dimensions
## do it properly (maybe via an auxiliary function)
wrk <- new("FlexibleLagged", data = pacr)
wrk[1:maxLag(wrk)]
}else {# lag_0 = "var" or anything else - put the variance in lag_0
## TODO: this may fail if 'x' doesn't contain info about autocovariances.
## Maybe just put NA instead?
wrk <- autocovariances(x, maxlag = 0, ...)
## BUGFIX: was: res <- as(obj, "FlexibleLagged")
res <- new("FlexibleLagged", data = pacr)
res[[0]] <- wrk[[0]]
res
}
}
setGeneric("partialAutocorrelations", signature = c("x", "maxlag", "lag_0") )
partialAutocovariances <- function(x, maxlag, ...){
# TODO: see handwritten notes "parcor"!
# TODO: for now don't use lag_0 since it is not supported consistently yet.
# wrk <- autocorrelations(x, lag_0 = "var", ...) # TODO: put R(0) at lag(0)
# TODO: only univariate for now
acvfobj <- if(missing(maxlag))
autocovariances(x, ...)
else
autocovariances(x, maxlag, ...)
pacvf <- .comboAcvf(acvfobj, "pacvf")
if(is(acvfobj, "SampleAutocovariances")){
new("SamplePartialAutocovariances", data = pacvf, as(acvfobj, "Fitted"))
}else
new("PartialAutocovariances", data = pacvf)
}
setGeneric("partialAutocovariances", signature = c("x", "maxlag") )
partialVariances <- function(x, ...){ # improvement? was: "predictionVariances"
# TODO: see handwritten notes "parcor"!
# see also comments in the default for partialAutocovariances()
# TODO: only univariate for now
acvfobj <- autocovariances(x, ...)
wrk <- .comboAcvf(acvfobj, "psigma2")
new("PartialVariances", data = wrk)
if(is(acvfobj, "SampleAutocovariances")){
new("SamplePartialVariances", data = wrk, as(acvfobj, "Fitted"))
}else
new("PartialVariances", data = wrk)
}
setGeneric("partialVariances") # was: setGeneric("predictionVariances")
backwardPartialVariances <- function(x,...) stop(x)
setGeneric("backwardPartialVariances") # setGeneric("backwardPredictionVariances")
partialCoefficients <-
backwardPartialCoefficients <- function(x, p) stop(x)
setGeneric("partialCoefficients") # setGeneric("predictionCoefficients")
setGeneric("backwardPartialCoefficients") # setGeneric("backwardPredictionCoefficients")
setMethod("autocovariances",
signature(x = "VirtualAutocovariances"),
function (x, ...)
{
stop("Can't compute autocovariances from objects from class ", class(x))
}
)
setMethod("autocovariances",
signature(x = "Autocovariances", maxlag = "missing"),
function (x, ...)
{
x
}
)
setMethod("autocovariances",
signature(x = "Autocovariances", maxlag = "ANY"),
function (x, maxlag, lag_0, ...)
{
x[] <- x[0:maxlag] # may introduce NA's
x
}
)
## TODO: set default value for maxlag?
setMethod("autocovariances",
signature(x = "VirtualArmaModel"),
function (x, maxlag, ...)
{
co <- modelCoef(x, convention = "BJ", ...)
sigma2 <- sigmaSq(x) # method for this how to call it:
# sigmaSq(), innovationVariance(), innovationVar()?
if(missing(maxlag))
maxlag <- max(length(co$ar), length(co$ma))
## TODO: compare to
## FitARMA::TacvfARMA(phi = co$ar, theta = co$theta, lag.max = maxlag)
res <- ltsa::tacvfARMA(phi = co$ar, theta = co$ma, maxLag = maxlag,
sigma2 = sigma2)
as(res, "Autocovariances")
}
)
setMethod("autocorrelations",
signature(x = "ANY", maxlag = "ANY", lag_0 = "ANY"),
function (x, maxlag, lag_0, ...)
{
obj <- autocorrelations(x, maxlag, ...)
if(isTRUE(lag_0))
obj
else if(identical(lag_0, FALSE)){
## drop lag 0 and the class
# if(!is.null(dim(obj@data)) && dim(obj@data) == 3)
# obj@data[ -1, , ]
# else
# obj@data[-1]
obj[seq(length = maxLag(obj))] # was: obj[][-1]
}else {# lag_0 = "var" or anything else - put the variance in lag_0
## TODO: this may fail if 'x' doesn't contain info about autocovariances.
## Maybe just put NA instead?
wrk <- autocovariances(x, maxlag = 0, ...)
res <- as(obj, "FlexibleLagged")
res[[0]] <- wrk[[0]]
res
}
}
)
setMethod("autocorrelations",
signature(x = "ANY", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
if(isS4(x))
stop("there is no applicable method for object from class ", class(x))
if(is.list(x) && all(names(x) %in% c("ar", "ma", "sigma2"))){
## theoretical ARMA autocov
wrk <- ARMAacf(x$ar, x$ma, lag.max = maxlag, ...) # `...' here may be only
# lag.max TODO: take care
# of lag_0
new("Autocorrelations", data = wrk)
}else{ # numeric, "ts", "mts" etc. ...
## sample autocorrelations
wrk <- if(missing(maxlag))
acf(x, plot = FALSE, type = "correlation", ...)
else
acf(x, lag.max = maxlag, plot = FALSE, type = "correlation", ...)
.basecfclass2S4(wrk)
}
## TODO: RESULT CHANGED FROM THAT IN pcts, ACCOMODATE FOR THIS IN pcts!
##
## res <- acfbase2sl(res) # convert to season/lag format
## res <- slMatrix(res)
## if(isTRUE(lag_0))
## res
## else if(identical(lag_0, FALSE))
## as.matrix(res)[-1]
## else {# lag_0 = "var" or anything else - put the variance in lag_0
## stop('lag_0 = "var" not implemented yet for this method.')
## }
}
)
setMethod("autocorrelations",
signature(x = "Autocorrelations", maxlag = "missing", lag_0 = "missing"),
function (x, ...)
{
## as(x, "Autocorrelations") # x may be from a subclass
x
}
)
setMethod("autocorrelations",
signature(x = "Autocorrelations", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, lag_0, ...)
{
## new("Autocorrelations", data = x[0:maxlag]) # may introduce NA's
x[] <- x[0:maxlag] # may introduce NA's
x
}
)
setMethod("autocorrelations",
signature(x = "Autocovariances", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
if(missing(maxlag))
maxlag <- maxLag(x)
acr <- x[0:maxlag] / x[[0]] # TODO: Only univariate case for now!
if(is(x, "Fitted"))
new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
else
new("Autocorrelations", data = acr)
}
)
## .pacvf2acrf <-
## function (x, maxlag, ...)
## {
## acr <- modelCoef(x, "Autocorrelations")
## res <- if(is(x, "Fitted"))
## new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
## else
## new("Autocorrelations", data = acr)
## if(!missing(maxlag))
## res@data <- res[0:maxlag] # may introduce NA's
## res
## }
##
## setMethod("autocorrelations",
## signature(x = "PartialAutocovariances", maxlag = "ANY", lag_0 = "missing"),
## .pacvf2acrf
## )
setMethod("autocorrelations",
signature(x = "PartialAutocorrelations", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
acr <- modelCoef(x, "Autocorrelations")
res <- if(is(x, "Fitted"))
new("SampleAutocorrelations", data = acr, as(x, "Fitted"))
else
new("Autocorrelations", data = acr)
if(!missing(maxlag))
## BugFix 2020-02-29 was: res@data <- res[0:maxlag]
## TODO: needs more testing
res@data <- Lagged(res[0:maxlag]) # may introduce NA's
res
}
)
## TODO: set default value for maxlag?
setMethod("autocorrelations",
signature(x = "VirtualArmaModel", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
## lazy - calls the default method with "BD" convention
co <- modelCoef(x, convention = "BD", ...)
co$sigma2 <- 1 # for autocorrelations this doesn't matter, but may be NA in x
if(missing(maxlag))
maxlag <- max(length(co$ar), length(co$ma))
autocorrelations(co, maxlag = maxlag, ...)
}
)
setMethod("autocorrelations",
signature(x = "VirtualSarimaModel", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
arma <- as(x, "ArmaModel") # should give error if x is non-stationary
autocorrelations(arma, maxlag = maxlag, ...)
}
)
setMethod("partialAutocorrelations",
signature(x = "ts", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...)
{
wrk <- if(missing(maxlag))
pacf(x, plot = FALSE)
else
pacf(x, lag.max = maxlag, plot = FALSE)
.basecfclass2S4(wrk)
}
)
setMethod("partialAutocorrelations",
signature(x = "mts", maxlag = "ANY", lag_0 = "missing"),
function (x, maxlag, ...){ # TODO: check what definition they use, etc.
wrk <- if(missing(maxlag))
pacf(x, plot = FALSE)
else
pacf(x, lag.max = maxlag, plot = FALSE)
r0 <- acf(x, lag.max = 0, plot = FALSE)
res <- .basecfclass2S4(wrk, r0)
## 2020-02-29 not needed, maxlag has already been used above!
## Also, this is wrong since res@data has a Lagged class.
## if(!missing(maxlag))
## res@data <- res[0:maxlag] # may introduce NA's
res
}
)
## BUGFIX: this cannot be done, see details elsewhere in this document!
## setMethod("partialAutocorrelations",
## signature(x = "PartialAutocovariances", maxlag = "ANY", lag_0 = "missing"),
## function (x, maxlag, ...)
## {
## pacr <- modelCoef(x, "PartialAutocorrelations")
## res <- new("PartialAutocorrelations", data = pacr)
## if(!missing(maxlag))
## res@data <- res[0:maxlag] # may introduce NA's
## res
## }
## )
## not needed since the default does the job
##
## setMethod("partialAutocorrelations",
## signature(x = "VirtualArmaModel"),
## function (x, maxlag, lag_0 = TRUE, ...)
## {
## co <- armaCoef(x, convention = "BJ", ...)
## sigma2 <- sigmaSq(x)
##
## wrk <- FitARMA::TacvfARMA(phi = numeric(0), theta = numeric(0), lag.max = maxlag)
## ## TODO: mult by sigma2
## }
## )
.comboAcvf <- function(acvf, ind){ # acvf is a Lagged object (actually Autocovariance or similar)
acvf <- acvf[] # drop the "Laggedness" (i.e. get the raw data)
## scalar case only; for now
r0 <- acvf[1]
acf <- acvf/r0
acfonly <- as(acf, "vector")[-1]
wrk <- ltsa::DLAcfToAR(acfonly) # further args: (useC = TRUE, PDSequenceTestQ = FALSE)
## wrk is a matrix with columns containing ar, pacf, sigma2;
## (but sigma2 (column "sigsqk") contains standardised innov vars,
## i.e. ones corresponding to R[0]=1);
## modify wrk below to our requirements.
## innovation variances in wrk are for R[0] = 1, so multiple them by R[0]
wrk[ , 3] <- r0 * wrk[ , 3]
wrk[ , 2] <- wrk[ , 2] * wrk[ , 3] # pacf => pacvf
res <- rbind(c(1, r0, r0), wrk) # prepend a row (1 R[0] R[0])'
# 2016-11-24 was: res <- cbind(acvf@data, res) # prepend also the acvf's
res <- cbind(acvf[], res) # prepend also the acvf's
# 'p' in 'pacvf' and 'psigma2' is for 'partial' as in pacf, etc.
colnames(res) <- c("acvf", "ar", "pacvf", "psigma2")
rownames(res) <- 0 : maxLag(acvf)
res <- t(res) # 2016-12-07 now return the lags as columns!
if(missing(ind))
res
else
res[ind, ]
}
## not used?
.comboAcvfNames <- c("Autocovariances" = "acvf",
"ar" = "ar", # not used
"PartialAutocovariances" = "pacvf",
"PartialVariances" = "sigma2" )
.comboAcf <- function(acvf, ind){ # acvf or acf
acvf <- acvf[] # drop the "Laggedness" (i.e. get the raw data)
## scalar case only; for now
r0 <- acvf[1]
acf <- acvf/r0
acfonly <- acf[-1]
wrk <- ltsa::DLAcfToAR(acfonly) # further args: (useC = TRUE, PDSequenceTestQ = FALSE)
## wrk is a matrix with columns containing ar, pacf, sigma2;
## (but sigma2 (column "sigsqk") contains standardised innov vars,
## i.e. ones corresponding to R[0]=1);
## modify wrk below to our requirements.
res <- cbind(acfonly, wrk) # prepend also the acf's
res <- rbind(c(1, 1, 1, 1), res) # prepend a row of 1's
colnames(res) <- c("acf", "ar", "pacf", "stdsigma2")
rownames(res) <- 0 : maxLag(acvf)
res <- t(res) # 2016-12-07 now return the lags as columns!
if(missing(ind))
res
else
res[ind, ]
}
## setAs("Autocovariances", "ComboAutocovariances",
## function(from){
##
##
## new("ComboAutocovariances", coef = )
## }
## )
setAs("vector", "Autocovariances", function(from) new("Autocovariances", data = from))
setAs("vector", "Autocorrelations", function(from) new("Autocorrelations", data = from))
setAs("vector", "PartialAutocovariances",
function(from) new("PartialAutocovariances", data = from))
setAs("vector", "PartialAutocorrelations",
function(from) new("PartialAutocorrelations", data = from))
## setAs("ANY", "PartialAutocovariances",
## function(from){
## partialAutocovariances(from)
## })
##
## setAs("ANY", "PartialVariances",
## function(from){
## partialVariances(from)
## })
##
## setAs("ANY", "Autocorrelations",
## function(from){
## autocorrelations(from)
## })
##
## setAs("ANY", "PartialAutocorrelations",
## function(from){
## partialAutocorrelations(from)
## })
setAs("ANY", "ComboAutocovariances",
function(from){
acvf <- as(from, "Autocovariances")
wrk <- .comboAcvf(acvf)
new("ComboAutocovariances", data = wrk)
})
setAs("ANY", "ComboAutocorrelations",
function(from){
acf <- as(from, "Autocorrelations")
wrk <- .comboAcf(acf)
new("ComboAutococorrelations", data = wrk)
})
setAs("Autocovariances", "ComboAutocovariances",
function(from){
wrk <- .comboAcvf(from)
new("ComboAutocovariances", data = wrk)
})
setAs("Autocovariances", "ComboAutocorrelations",
function(from){
wrk <- .comboAcf(from)
new("ComboAutocorrelations", data = wrk)
})
setAs("Autocorrelations", "ComboAutocovariances",
function(from){
stop("not possible, missing R(0)")
})
setAs("Autocorrelations", "ComboAutocorrelations",
function(from){
wrk <- .comboAcf(from)
new("ComboAutocorrelations", data = wrk)
})
## setAs("PartialVariances", "Autocovariances",
## function(from){
## stop("PartialVariances cannot be converted to the requested class.")
## })
##
## setAs("PartialVariances", "ComboAutocovariances",
## function(from){
## stop("PartialVariances cannot be converted to the requested class.")
## })
##
## setAs("PartialVariances", "Autocorrelations",
## function(from){
## stop("PartialVariances cannot be converted to the requested class.")
## })
##
## setAs("PartialVariances", "ComboAutocorrelations",
## function(from){
## stop("PartialVariances cannot be converted to the requested class.")
## })
setAs("Autocorrelations", "Autocovariances",
function(from){
new("Autocovariances", data = from[])
})
setAs("SampleAutocorrelations", "SampleAutocovariances",
function(from){
new("SampleAutocovariances", data = from[], as(from, "Fitted"))
})
setMethod("modelCoef", c("VirtualAutocovariances", "missing", "missing"),
function(object){
object@data
}
)
setMethod("modelCoef", c("VirtualAutocovariances", "character", "missing"),
function(object, convention){
if(identical(class(object), convention))
## this case is equivalent to the one-argument call()
modelCoef(object)
else{
convention <- new(convention)
modelCoef(object, convention = convention )
}
}
)
setMethod("modelCoef", c("VirtualAutocovariances", "VirtualAutocovariances", "missing"),
function(object, convention){
acvf <- modelCoef(object, "Autocovariances")
modelCoef(acvf, convention)
}
)
setMethod("modelCoef", c("VirtualAutocovariances", "Autocovariances", "missing"),
function(object, convention){
stop("Can't obtain autocovariances from object from class ", class(object))
}
)
setMethod("modelCoef", c("Autocovariances", "ComboAutocovariances", "missing"),
function(object, convention){
.comboAcvf(object)
}
)
setMethod("modelCoef", c("Autocovariances", "ComboAutocorrelations", "missing"),
function(object, convention){
.comboAcf(object)
}
)
# not defined since R[0] is unknown
# setMethod("modelCoef", c("Autocorrelations", "ComboAutocovariances", "missing") )
setMethod("modelCoef", c("Autocorrelations", "ComboAutocorrelations", "missing"),
function(object, convention){
.comboAcf(object)
}
)
setMethod("modelCoef", c("ComboAutocovariances", "Autocovariances", "missing"),
function(object, convention){
## 2020-02-29 changing this (here and at similar places) with temporary patch
## since Lagged2d objects don't have method for "character" (also "logical")
##
## TODO: UPDATE PACKAGE "lagged" with such methods, if sensible!
##
## object@data["acvf", ]
tmp <- which(rownames(object@data[]) == "acvf")
object@data[tmp, ]
}
)
## TODO: to AR model (or filter) using $ar
setMethod("modelCoef", c("ComboAutocovariances", "PartialAutocovariances", "missing"),
function(object, convention){
## object@data["pacvf", ]
tmp <- which(rownames(object@data[]) == "pacvf")
object@data[tmp, ]
}
)
setMethod("modelCoef", c("ComboAutocovariances", "PartialVariances", "missing"),
function(object, convention){
## object@data["sigma2", ]
tmp <- which(rownames(object@data[]) == "sigma2")
object@data[tmp, ]
}
)
## acvf to acf etc.
setMethod("modelCoef", c("ComboAutocovariances", "VirtualAutocovariances", "missing"),
function(object, convention){
acvf <- modelCoef(object, "Autocovariances")
modelCoef(acvf, convention)
}
)
setMethod("modelCoef", c("ComboAutocorrelations", "Autocorrelations", "missing"),
function(object, convention){
## object@data["acf", ]
tmp <- which(rownames(object@data[]) == "acf")
object@data[tmp, ]
}
)
## TODO: to AR model (or filter) using $ar
setMethod("modelCoef", c("ComboAutocorrelations", "PartialAutocorrelations", "missing"),
function(object, convention){
## object@data["pacf", ]
tmp <- which(rownames(object@data[]) == "pacf")
object@data[tmp, ]
}
)
## need class StandardizedPartialVariances
## setMethod("modelCoef", c("ComboAutocorrelations", "StandardizedPartialVariances", "missing"),
## function(object, convention){
## object@data["stdsigma2", ]
## }
## )
setMethod("modelCoef", c("Autocovariances", "Autocorrelations", "missing"),
function(object, convention){
co <- object[] # or modelCoef(object) - no need for such generality?
co / co[1] # TODO: modify to work in multivariate case
}
)
setMethod("modelCoef", c("Autocovariances", "PartialAutocorrelations", "missing"),
function(object, convention){
# alternatively, obtain them from combo (which uses ltsa::DLAcfToAR;
# using here acf2AR; compare to combo during testing.
co <- object[] # or modelCoef(object) - no need for such generality?
ar.all <- acf2AR(co) # all partial prediction coef's (YW)
c(1, diag(ar.all) )
}
)
## same as for "Autocovariances"
setMethod("modelCoef", c("Autocorrelations", "PartialAutocorrelations", "missing"),
function(object, convention){
co <- object@data # or modelCoef(object) - no need for such generality?
ar.all <- acf2AR(co) # all partial prediction coef's (YW)
c(1, diag(ar.all) )
}
)
setMethod("modelCoef", c("PartialAutocorrelations", "Autocorrelations", "missing"),
function(object, convention){
# 2016-11-24 was:
# pacf <- object@data
# pacfonly <- pacf[-1]
pacfonly <- object[1:maxLag(object)]
# 2022-02-14 was: ar <- FitAR::PacfToAR(pacfonly)
ar <- pacf2Ar(pacfonly)
# alternatively:
# acvf <- ltsa::tacvfARMA(phi = ar, maxLag = length(ar))
# acf <- acvf/acvf[0]
ARMAacf(ar, lag.max = length(ar))
}
)
##not available on purpose
##setMethod("modelCoef", c("Autocorrelations", "PartialAutocovariances", "missing"))
## setMethod("modelCoef", c("PartialAutocovariances", "PartialAutocorrelations", "missing"),
## function(object, convention){
## ## pacr <- object[] / object[0] # TODO: currently scalar only
## stop("Partial autocorrelations are not uniquely determined by partial autocovariances")
## }
## )
##
## setMethod("modelCoef", c("PartialAutocovariances", "Autocorrelations", "missing"),
## function(object, convention){
## ## pacfonly <- object[] / object[0] # TODO: currently scalar only
## ## ar <- FitAR::PacfToAR(pacfonly[-1]) # drop lag 0 !
## ## # alternatively:
## ## # acvf <- ltsa::tacvfARMA(phi = ar, maxLag = length(ar))
## ## # acf <- acvf/acvf[0]
## ## ARMAacf(ar, lag.max = length(ar))
## stop("Autocorrelations are not uniquely determined by partial autocovariances")
## }
## )
## n - length of time series, scalar
## npar - no. estimated param, scalar
## nlags - number of acf for the LB test, vector of positive integers
##
acfIidTest <- function(acf, n, npar = 0, nlags = npar + 1,
method = c("LiMcLeod", "LjungBox", "BoxPierce"),
interval = 0.95, expandCI = TRUE, ...){
method <- match.arg(method)
maxlag <- max(nlags)
usedLags <- 1:maxlag
switch(method,
"LjungBox" = {
rsq <- acf[usedLags]^2 / (n - usedLags)
Q <- n * (n+2) * cumsum(rsq)
},
"BoxPierce" = {
rsq <- acf[usedLags]^2
Q <- n * cumsum(rsq)
},
"LiMcLeod" = {
rsq <- acf[usedLags]^2
Q <- n * cumsum(rsq) + usedLags * (usedLags + 1) / (2*n)
}#, stop("Unknown method")
)
Q <- Q[nlags] # keep only the requested ChiSq values
df <- nlags - npar
df[df <= 0] <- 1 # TODO: this should really be NA or depend on the estimation method
pval <- pchisq(Q, df = df, lower.tail = FALSE)
wrk <- cbind(ChiSq = Q, DF = df, pvalue = pval)
attr(wrk, "method") <- method
res <- list(test = wrk)
if(!is.null(interval)){
lq <- qnorm((1 - interval) / 2)
int <- lq / sqrt(n)
if(expandCI)
int <- rep(int, maxlag)
int <- cbind(int, - int)
res$ci <- int
attr(res$ci, "level") <- interval
}
res
}
setGeneric("acfIidTest")
setMethod("acfIidTest", "SampleAutocorrelations",
function(acf, n, npar, nlags, method, interval = 0.95, x){
if(missing(n))
n <- acf@n
else
if(n != acf@n){
## only warning to allow 'what-if use.
warning("argument 'n' is not equal to slot 'n' of the autocorrelations")
}
callNextMethod(acf, n = n, npar = npar, nlags = nlags, method = method, interval = interval)
})
setMethod("acfIidTest", "missing",
function(acf, n, npar, nlags, method, interval = 0.95, x){
if(missing(n))
n <- length(x)
acf <- autocorrelations(x, maxlag = max(nlags))
acfIidTest(acf, n = n, npar = npar, nlags = nlags, method = method, interval = interval)
})
acfMaTest <- function(acf, ma, n, nlags,
interval = 0.95){
if(!is(acf, "Lagged"))
acf <- Lagged(acf)
maxlag <- max(nlags)
usedLags <- 1:maxlag
W <- nvcovOfAcfBD(acf, ma = ma, maxlag = maxlag)
Q <- numeric(length(nlags))
for(i in seq(along = nlags)){
h <- nlags[i]
la <- (ma + 1):h
rho <- acf[la]
stat <- rho %*% solve(W[la,la]) %*% rho
Q[i] <- stat / n
}
df <- nlags - ma
pval <- pchisq(Q, df = df, lower.tail = FALSE)
wrk <- cbind(ChiSq = Q, DF = df, pvalue = pval)
attr(wrk, "method") <- paste0("Ma(", ma, ") Test")
res <- list(test = wrk)
if(!is.null(interval)){
lq <- qnorm((1 - interval) / 2, sd = sqrt(diag(W)))
int <- lq / sqrt(n)
## 2019-03-15 <BUGFIX> was: int <- cbind(int, - int)
## also need levels for the lags!
rho_hat_null <- c(acf[seq_len(ma)], rep(0, length(int) - ma))
int <- rho_hat_null + cbind(int, - int)
res$ci <- int
attr(res$ci, "level") <- interval
}
res
}
## 'x' is the time series
.wnacf_asycov_iid <- function(x, maxlag, forcor = TRUE, n = length(x)){
if(forcor)
diag(1/n, nrow = maxlag)
else{
r0 <- acf(x, type = "covariance", lag.max = 0, plot = FALSE)$acf[1, , ]
diag(r0^2/length(x), nrow = maxlag) # TODO: verify!
}
}
.wnacf_asycov_garch <- function(x, maxlag, forcor = TRUE){ # TODO: default for maxlag
x <- as.vector(x) # otherwise x*wrk below doesn't work
n <- length(x)
wrk <- sapply(1:maxlag, function(i) c(rep(0, i), x[1:(n-i)]))
wrk <- x * wrk # multiple each column by x using recycling
res <- acf(wrk, type = "covariance", lag.max = 0, plot = FALSE, demean = FALSE)
res <- res$acf[1, , ] # extract the zero lag matrix
if(forcor){
r0 <- acf(x, type = "covariance", lag.max = 0, plot = FALSE,
demean = FALSE)$acf[1, , ]
res <- res / r0^2
}
res / n # since sqrt(n)* rhat has asy.cov. 'res'
}
.wnacf_asyse <- function(x, maxlag, h0 = "iid", forcor = TRUE, n = length(x)){
switch(h0,
"iid" = {
v <- .wnacf_asycov_iid(x = x, maxlag = maxlag, forcor = forcor, n = n)
},
"garch" = {
v <- .wnacf_asycov_garch(x = x, maxlag = maxlag, forcor = forcor)
},
##default
stop("Unrecognised null hypothesis")
)
ses <- sqrt(diag(v))
lags <- 1:maxlag
cbind(Lag = lags, StdError = ses) # todo: include 'Estimate = cf[parm]'
}
acfGarchTest <- function(acr, x, nlags, interval = 0.95){
v <- .wnacf_asycov_garch(x, maxlag = max(nlags))
# the formulas in the book multiply by n, since Gamma is 'n times asy.cov.'
# n <- length(x)
# Q <- n * sapply(nlags, function(h) acr[1:h] %*% solve(v[1:h,1:h]) %*% acr[1:h] )
Q <- sapply(nlags, function(h) acr[1:h] %*% solve(v[1:h,1:h]) %*% acr[1:h] )
pval <- pchisq(Q, df = nlags, lower.tail = FALSE)
res <- list(test = cbind(h = nlags, Q = Q, pval = pval))
if(!is.null(interval)){
lq <- qnorm((1 - interval) / 2)
se <- sqrt(diag(v))
int <- lq * se
int <- cbind(int, -int)
res$ci <- int
}
res
}
acfWnTest <- function(acr, x, nlags, interval = 0.95, ...){
maxlag <- max(nlags)
n <- length(x)
v <- nvarOfAcfKP(x, maxlag = maxlag, ...)
## nvcov is diagonal in this case
Q <- n * cumsum(acr[1:maxlag]^2 / v )[nlags]
pval <- pchisq(Q, df = nlags, lower.tail = FALSE)
res <- list(test = cbind(h = nlags, Q = Q, pval = pval))
if(!is.null(interval)){
lq <- qnorm((1 - interval) / 2)
# se <- sqrt(diag(v))
se <- sqrt(v / n)
int <- lq * se
int <- cbind(int, -int)
res$ci <- int
}
res
}
whiteNoiseTest <- function(object, h0, ...){
switch(h0,
"iid" = {
acfIidTest(object, ...)
},
"garch" = {
acfGarchTest(object, ...) # no method yet for this case
},
"sv" = ,
"arch-type" = {
acfWnTest(object, ...)
},
##default
stop("Unrecognised null hypothesis")
)
}
setMethod("coef", "SampleAutocorrelations",
function (object, ...){
dataWithLagNames(object)
}
)
setMethod("vcov", "SampleAutocorrelations",
function (object, assuming = "iid", maxlag = maxLag(object), ...) {
# res <- switch(assuming,
# iid = .wnacf_asycov_iid(n = object@n, maxlag = maxlag),
# garch = .wnacf_asycov_garch(maxlag = maxlag, ...),
# stop("Unknown assumption")
# )
# if(is.null(names(res)))
# ## TODO: currently scalar case only;
# colnames(res) <- rownames(res) <- paste0("Lag_", 1:nrow(res))
# res
if(inherits(assuming, "Arima")){
## todo: could get vcov of params and transform it for acf's
## maybe have additional argument here (method?)
mo <- list(ar = assuming$model$phi,
ma = assuming$model$theta,
sigma2 = assuming$sigma2 )
res <- nvcovOfAcf(mo, maxlag = maxlag) / object@n
res
}else
vcovAcf(object, assuming, maxlag = maxlag, ...)
}
)
setGeneric("vcovAcf", function(object, assuming, ...) stop("no suitable method"))
setMethod("vcovAcf", c("SampleAutocorrelations", "character"),
function (object, assuming, maxlag = maxLag(object), ...) {
if(length(assuming) > 1)
stop("'assuming' must be of length 1")
res <- switch(assuming,
iid = .wnacf_asycov_iid(n = object@n, maxlag = maxlag),
garch = .wnacf_asycov_garch(maxlag = maxlag, ...),
stop("Unknown assumption")
)
if(is.null(names(res)))
## TODO: currently scalar case only;
colnames(res) <- rownames(res) <- paste0("Lag_", 1:nrow(res))
res
}
)
setMethod("vcovAcf", c("SampleAutocorrelations", "ArmaModel"),
function (object, assuming, maxlag = maxLag(object), ...) {
res <- nvcovOfAcf(c(modelCoef(assuming), list(sigma2 = sigmaSq(assuming))),
maxlag = maxlag) / object@n
res
}
)
se <- function (object, ...) {
vc <- vcov(object, ...)
sqrt(diag(vc))
}
## TODO: make S4 generic
.fixed_values_under_h0 <- function(object, lags){
m <- length(lags)
if(identical(object, "iid"))
numeric(m)
else if(identical(object, "garch"))
numeric(m)
else if(inherits(object, "Arima")){
order <- object$arma
names(order) <- c("ar", "ma", "sar", "sma", "period", "d", "ds")
if(order["d"] != 0 || order["ds"] != 0)
stop("Can't compute acf confint for non-stationary models")
else if(all(order[c("ar", "sar")] == 0)){
if(order["sma"] == 0){
## TODO: can be refined
res <- ifelse(lags <= order["ma"], NA_real_, 0)
res
}else if(order["ma"] == 0){ # pure seasonal MA
res <- ifelse(lags %% order["period"] == 0 &
lags %/% order["period"] <= order["ma"] ,
NA_real_, 0)
res
}else
res <- rep(NA_real_, m)
}
}else if(is(object, "ArmaModel")){
## for a theoretical model all acfs are fixed
## for fitted models => see the case 'Arima'
## TODO: need further methods here
# if(modelOrder(object)$ar == 0)
# ifelse(lags <= modelOrder(object)$ma, NA_real_, 0)
# else
#browser()
autocovariances(c(modelCoef(object), list(sigma2 = sigmaSq(object))),
maxlag = max(lags))[lags]
}else ## default method - all not fixed
rep(NA_real_, m)
}
# setMethod("diagOfVcov", "SampleAutocorrelations",
setMethod("confint", "SampleAutocorrelations",
function (object, parm, level = 0.95, se = FALSE, maxlag, ..., assuming){
## based on confint.default(), main difference: we pass "..." to vcov();
## arg. 'se'
cf <- coef(object)[-1] # drop lag 0
pnames <- names(cf)
if (missing(parm)){
if(missing(maxlag)){
lags <- 1:length(pnames)
parm <- pnames
}else{
lags <- 1:maxlag
parm <- pnames[1:maxlag]
}
}else if (is.numeric(parm)){
lags <- parm
parm <- pnames[parm]
}
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- .stats.format.perc(a, 3)
fac <- qnorm(a)
ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm, pct))
# TODO: pass maxlag to vcov?
# 2022-01-17 was:
# v <- vcov(object, ...)
# ses <- sqrt(diag(v))[parm]
ses <- se(object, assuming = assuming, ...)[lags]
cg <- cf[parm]
if(!is.null(assuming)){
int_types <- .fixed_values_under_h0(assuming, lags)
if(!all(is.na(int_types))){
cg[!is.na(int_types)] <- int_types[!is.na(int_types)]
}
}else
int_types <- rep(NA_real_, length(parm))
ci[] <- cg + ses %o% fac
ci <- cbind(Lag = lags, ci, Estimate = cf[parm], H0_fixed = int_types)
if(se)
ci <- cbind(ci, StdError = ses)
attr(ci, "assuming") <- assuming
ci
}
)
.plot.acf.test <- function (x, y, data, method = "LiMcLeod", nlags, interval = 0.95,
legend = TRUE, ylim = c(NA, NA), ylim.fac = 1.2,
xlab = "Lag",
ci.lty = 2:3, # 2017-08-29 new, corresponding changes below
h0 = "garch", # 2019-04-09 new to accommodate the new white noise test
# (Kokoszka-Politis)
...){
if(missing(nlags))
nlags <- seq(from = maxLag(x), to = 1, by = -5)
x.iid <- whiteNoiseTest(x, h0 = "iid", n = x@n, interval = interval, # x = data,
method = method, nlags = nlags)
ci <- x.iid$ci
y <- cbind(1:nrow(ci), x.iid$ci, Estimate = x[1:nrow(ci)])
ylim.loc <- ylim
if(is.na(ylim[1])) ylim.loc[1] <- min(y[, 2])
if(is.na(ylim[2])) ylim.loc[2] <- max(y[, 3])
if(!missing(data)){
x.garch <- whiteNoiseTest(x, h0 = "garch", x = data, nlags = nlags,
interval = interval)
ylim <- c(min(ylim[1], x.garch$ci[ , 1]),
max(ylim[2], x.garch$ci[ , 2])
)
## modify 'ylim.loc' if the corresponding 'ylim' element is NA
if(is.na(ylim[1])) ylim.loc[1] <- min(ylim.loc[1], x.garch$ci[ , 1])
if(is.na(ylim[2])) ylim.loc[2] <- max(ylim.loc[2], x.garch$ci[ , 2])
}else{ # todo: modify some parameters for this case.
legend <- FALSE
}
## replace NA's in 'ylim' with the computed values
if(is.na(ylim[1])) ylim[1] <- ylim.fac * ylim.loc[1]
if(is.na(ylim[2])) ylim[2] <- ylim.fac * ylim.loc[2]
.ciplot(x = y,
## lag_0 = TRUE,
ylim = ylim,
ci.col = "brown",
ci.lty = ci.lty[1],
xlab = xlab, ## xlab = colnames(x)[1],
## ylab = "Estimate & CI",
## ann = TRUE,
...
)
if(!missing(data)){
x.garch <- whiteNoiseTest(x, h0 = h0, x = data, nlags = nlags)
.cilines(1:nrow(x.garch$ci), x.garch$ci, col = "blue", lty = ci.lty[2])
}
## see ?legend, example thanks to Uwe Ligges
if(legend){
legend("topright", legend = c("H0: iid", paste0("H0: ", h0)),
col = c("brown ", "blue"), lty = ci.lty )
}
invisible()
}
.cilines <- function(x, y, type = "l", col = "blue", lty = 2, ...){
lines(x, y[ , 1], type = type, col = col, lty = lty, ...)
lines(x, y[ , 2], type = type, col = col, lty = lty, ...)
invisible()
}
.plot_acr <- function(x, y, xlim = c(0,max(x)), ylim = c(-1, 1), type = "h", ...){
plot(x, y, xlim = xlim, ylim = ylim, type = type, # xlab = xlab, ylab = ylab,
...)
abline(h=0)
## if(ann)
## title(xlab = xlab, ylab = ylab, ...)
## else
## title(...)
invisible()
}
.ciplot <- function(x, y, lag_0 = TRUE, ylim = c(-1, 1),
ci.col = "blue", ci.lty = 2,
xlab = colnames(x)[1],
ylab = "Estimate & rejection levels",
ann = TRUE,
...){
lags <- xval <- x[ , 1]
yval <- x[ , "Estimate"]
## lbval <- x[ , 2]
## ubval <- x[ , 3]
## TODO: need more care here.
if(lag_0){
xval <- c(0, xval)
yval <- c(1, yval)
}
## plot(numeric(0), xlim = c(0,max(xval)), ylim = ylim, xlab = xlab, ylab = ylab,
## ann = FALSE,
## ...)
## abline(h=0)
## lines(c(0,0), c(0,1))
.plot_acr(xval, yval, xlim = c(0,max(xval)), ylim = ylim, xlab = xlab, ylab = ylab,
ann = ann, ...)
## lines(xval, yval, type = "h")
.cilines(lags, x[ , 2:3], type = "l", col = ci.col, lty = ci.lty)
## if(ann)
## title(xlab = xlab, ylab = ylab, ...)
## else
## title(...)
invisible()
}
setMethod("plot", c(x = "SampleAutocorrelations", y = "matrix"),
function (x, y, main = "Acf test", ...){
.ciplot(x = y,
lag_0 = FALSE, ## ylim = c(-1, 1),
main = main,
## ci.col = "blue", ci.lty = 2,
## xlab = colnames(x)[1],
## ylab = "Estimate & CI",
## ann = TRUE,
...)
invisible()
}
)
setMethod("plot", c(x = "SampleAutocorrelations", y = "missing"),
function (x, y,
# data, method = "LiMcLeod", nlags, legend = TRUE,
main = "Acf test",
...){
.plot.acf.test(x = x, main = main, ...)
}
)
setMethod("plot", c(x = "SamplePartialAutocorrelations", y = "missing"),
function (x, y,
# data, method = "LiMcLeod", nlags, legend = TRUE,
main = "Pacf test",
...){
.plot.acf.test(x = x, main = main, ...)
}
)
acfOfSquaredArmaModel <- function(model, maxlag){
ar <- if(is.null(model$ar)) numeric(0) else model$ar
arsq <- - coef(polynom(c(1, - ar))^2)[-1]
masq <- coef(polynom(c(1, model$ma))^2)[-1]
Sg <- model$sigma2 ^ 2
autocovariances(list(ar = arsq, ma = masq, sigma2 = Sg), maxlag = maxlag)
}
## naive implementation
nvcovOfAcf <- function(model, maxlag){
R <- autocovariances(model, maxlag = maxlag)
r <- R / R[0]
Rg <- acfOfSquaredArmaModel(model, maxlag = 2 * maxlag)
res <- matrix(NA_real_, nrow = maxlag, ncol = maxlag)
for(k in 1:maxlag){
for(l in k:maxlag){
res[k, l] <- Rg[l - k] + Rg[l + k] - 2 * Rg[k] * r[l] -
2 * Rg[l] * r[k] + 2 * Rg[0] * r[k] * r[l]
if(k != l)
res[l, k] <- res[k, l]
}
}
res / R[0]^2
}
nvcovOfAcfBD <- function(acf, ma = NULL, maxlag){
maacf <- new("Lagged1d", data = acf[])
max_avail_lag <- maxLag(maacf)
if(is.null(ma))
ma <- max_avail_lag
else
stopifnot(ma <= max_avail_lag)
## TODO: this may make max_avail_lag too large. Fix!
if(max_avail_lag < 2 * maxlag + ma)
max_avail_lag <- 2 * maxlag + ma
if(ma < max_avail_lag)
maacf[(ma + 1) : max_avail_lag] <- 0
## eq. 7.2.6. BD, p. 222
res <- matrix(NA_real_, nrow = maxlag, ncol = maxlag)
for(j in seq_len(maxlag)){
for(i in j:maxlag){ # i >= j
maxk <- ma + max(i, ma)
val <- 0
for(k in 1 : maxk){
## k + i can be up to ma + 2*maxlag
val <- val + (maacf[k + i] + maacf[abs(k - i)] - 2 * maacf[i] * maacf[k]) *
(maacf[k + j] + maacf[abs(k - j)] - 2 * maacf[j] * maacf[k])
}
res[i, j] <- val
if(i != j)
res[j, i] <- res[i, j]
}
}
res
}
nvarOfAcfKP <- function(x, maxlag, center = FALSE, acfscale = c("one", "mom")){
acfscale <- match.arg(acfscale)
if(center)
x <- x - mean(x)
acv <- acf(x^2, lag.max = maxlag, demean = FALSE, type = "covariance", plot = FALSE)
acv <- acv$acf[ , , ] # make numeric, or more explicitly: acv$acf[ , , 1]
## the denominator, NOTE: it is the square of the lag 0 acv of x, not x^2
n <- length(x)
den <- (sum(x^2) / n) ^ 2
switch(acfscale,
one = acv[-1] / den,
mom = {
fac <- n / (n - seq_len(maxlag))
fac * acv[-1] / den
}
## 2020-02-29 removing this since match.arg() above catches this error.
## ,
## ## default
## stop("argument 'acfscale' must be one of 'one', 'mom' or abbreviation thereof")
)
}
rgarch1p1 <- function(n, alpha, beta, omega, n.skip = 100){
Esigma2 <- omega / (1 - alpha - beta)
ntot <- n.skip + n
eta <- rnorm(ntot)
sigma2 <- eps <- numeric(ntot)
sigma2[1] <- Esigma2 # use E sigma_t^2 as initial value
eps[1] <- sqrt(sigma2[1]) * eta[1]
for(t in 2:ntot){
sigma2[t] <- omega + alpha * eps[t-1]^2 + beta * sigma2[t-1]
eps[t] <- sqrt(sigma2[t]) * eta[t]
}
eps[(n.skip + 1):ntot]
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.