Nothing
#Source: artfimaTACVF.R
#autocovariance for stationary ARTFIMA
#
"artfimaTACVF" <-
function(d=numeric(0), lambda=numeric(0), phi = numeric(0),
theta = numeric(0), maxlag, sigma2 = 1, obj=NULL)
{
if (!is.null(obj)) {
if ("artfima" == class(obj)) {
d <- obj$dHat
lambda <- obj$lambdaHat
phi <- obj$phiHat
theta <- obj$thetaHat
sigma2 <- obj$sigmaSq
}
}
#**WARNING**:assumes phi and theta in admissible region!!
#use invertibleQ(phi) and invertibleQ(theta) prior to using artfimaTACVF!
ARMALength <- sum(length(phi),length(theta))
ARTFIMALength <- ARMALength+sum(length(d), length(lambda))
isWhiteNoise <- identical(0==sum(ARTFIMALength),TRUE)
if (isWhiteNoise) return(c(sigma2, rep(0, maxlag)))
#pure ARMA
isARMA <- identical(length(d)==0,TRUE)||identical(d==0,TRUE)
if (isARMA)
return(tacvfARMA(phi = phi, theta = theta, maxlag = maxlag,
sigma2 = sigma2))
#fractional acvf
lagTrunc <- 2*max(128, nextn(maxlag, factors = 2))
if (identical(lambda>1e-7,TRUE)) {
x <- tacvfFI(d = d, lambda=lambda, maxlag = lagTrunc)
} else {
x <- tacvfFDWN(dfrac = d, maxlag = lagTrunc)
}
if (ARMALength==0) return(sigma2*x[1:(maxlag+1)])
#ARMA case
y <- tacvfARMA(phi = phi, theta = theta, maxlag = lagTrunc, sigma2 = 1)
z <- sigma2*mix(x, y)
z[1:(maxlag+1)]
}
#
#note domain of d are the real numbers
#requires package gsl
tacvfFI <-
function(d, lambda, maxlag, sigma2=1){
if (abs(d)<1e-8) return(c(1,rep(0,maxlag)))
k <- 0:maxlag
if (d>0) {
exL <- min(exp(-2*lambda), 0.99)
A <- hyperg_2F1(d, d + k, 1 + k, exL)
#HACK - handle error with hyper_2F1() !!!!!!!!!!!!!!!!!!!!
if (NaN %in% A) {#check for error
stop("error: NaN from hyperg_2F1()in tacvfFI()")
}
C <- k*lambda + lgamma(1+k)
B <- lnpoch(d,k)
ans <- A*exp(B-C)
} else { #Approximation when d<0. When d<-0.499, acvf is neglible.
ans <- exp(-lambda*(0:maxlag))*tacvfFDWN(max(d,-0.499), maxlag)
}
sigma2*ans
}
#
"tacvfFARMA" <-
function(phi = numeric(0), theta = numeric(0), dfrac = 0,
maxlag, sigma2 = 1)
{
#phi, theta: NOT VALIDATED!!
#assumes phi and theta in admissible region and abs(d)<0.5
#white noise case
ARMALength <- sum(length(phi),length(theta))
isWhiteNoise <- 0==sum(ARMALength,dfrac) #works if drfrac=0 or numeric(0)
if (isWhiteNoise) return(c(sigma2, rep(0, maxlag)))
isARMA <- identical(0==dfrac,TRUE)||(length(dfrac)==0)
if (isARMA) return(tacvfARMA(phi = phi, theta = theta, maxlag = maxlag,
sigma2 = sigma2))
lagTrunc <- 2*max(128, nextn(maxlag, factors = 2))
x <- tacvfFDWN(dfrac = dfrac, maxlag = lagTrunc)
if (ARMALength==0) return(sigma2*x[1:(maxlag+1)])
y <- tacvfARMA(phi = phi, theta = theta, maxlag = lagTrunc, sigma2 = 1)
z <- sigma2*mix(x, y)
z[1:(maxlag+1)]
}
#acvf fractionally differenced white noise
#
"tacvfFDWN" <-
function(dfrac, maxlag, sigma2=1)
{
if (dfrac>0.499) dfrac <- 0.499
x <- numeric(maxlag + 1)
x[1] <- gamma(1 - 2 * dfrac)/gamma(1 - dfrac)^2
for(i in 1:maxlag)
x[i + 1] <- ((i - 1 + dfrac)/(i - dfrac)) * x[i]
x*sigma2
}
"tacvfARMA" <-
function(phi = numeric(0), theta = numeric(0), maxlag = 20,
sigma2 = 1)
{
p <-length(phi)
q <- length(theta)
maxlagp1 <- maxlag + 1
if(max(p, q) == 0) {
return(c(sigma2, numeric(maxlag)))
}
r <- max(p, q) + 1
b <- numeric(r)
C <- numeric(q + 1)
C[1] <- 1
theta2 <- c(-1, theta)
phi2 <- numeric(3 * r)
phi2[r] <- -1
if(p > 0) {
phi2[r + 1:p] <- phi
}
if(q > 0) {
for(k in 1:q) {
C[k + 1] <- (-theta[k])
if(p > 0) {
for(i in 1:min(p, k)) {
C[k + 1] <- C[k + 1] + phi[i] * C[k + 1 - i]
}
}
}
}
for(k in 0:q) {
for(i in k:q) {
b[k + 1] <- b[k + 1] - theta2[i + 1] * C[i - k + 1]
}
}
if(p == 0) {
g <- c(b, numeric(maxlagp1))[1:maxlagp1]
return(g)
}
else if(p > 0) {
a <- matrix(numeric(r^2), ncol = r)
for(i in 1:r) {
for(j in 1:r) {
if(j == 1) {
a[i, j] <- phi2[r + i - 1]
}
else if(j != 1) {
a[i, j] <- phi2[r + i - j] + phi2[r + i + j - 2]
}
}
}
g <- solve(a, - b)
if(length(g) <= maxlag) {
g <- c(g, numeric(maxlagp1 - r))
for(i in (r + 1):maxlagp1) {
g[i] <- phi %*% g[i - 1:p]
}
return(sigma2*g[1:maxlagp1])
}
else if(length(g) >= maxlagp1) {
return(sigma2*g[1:maxlagp1])
}
}
}
"symtacvf" <-
function(x)
{
c(rev(x[-1])[-1], x)
}
mix <- function(x, y) {
n <- 2*length(x)-2
rev(Re(
fft(fft(symtacvf(x)) * fft(symtacvf(y)), inverse = TRUE)/n
)[(n/2 - 1):(n - 1)])
}
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.