R/tacvfARMA.R In ltsa: Linear Time Series Analysis

````tacvfARMA` <-
function(phi = numeric(0), theta = numeric(0), maxLag = 1, sigma2=1)
{
stopifnot(maxLag >= 0, sigma2>0)
#check for stationarity
p <- length(phi)
if(p>0){
pi=numeric(p)
phik <- phi
for (k in 1:p){
L <- p+1-k
a <- phik[L]
pi[L+1-k] <- a
phikp1 <- phik[-L]
if(is.na(a) || abs(a)==1)
stop("error: roots outside unit circle -- nonstationary/noncausal model")
phik <- (phikp1+a*rev(phikp1))/(1-a^2)
}
if (!all(abs(pi)<1))
stop("error: roots outside unit circle -- nonstationary/noncausal model")
}
#model is stationary, compute acvf
q <- length(theta)
if(max(p, q) == 0)
return(c(sigma2, numeric(maxLagp1))) #white noise case
maxLagp1 <- maxLag + 1
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]
else {
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
a[i, j] <- phi2[r + i - j] + phi2[r + i + j - 2]
}
g <- solve(a,  - b)
if(length(g) <= maxLag) {
g <- c(g, numeric(maxLag - r))
for(i in (r + 1):maxLagp1)
g[i] <- phi %*% g[i - 1:p]
}
}
sigma2*g[1:maxLagp1]
}
```

Try the ltsa package in your browser

Any scripts or data that you put into this service are public.

ltsa documentation built on May 2, 2019, 4:01 a.m.