| # File src/library/stats/R/ARMAtheory.R |
| # Part of the R package, https://www.R-project.org |
| # |
| # Copyright (C) 1995-2015 The R Core Team |
| # |
| # This program is free software; you can redistribute it and/or modify |
| # it under the terms of the GNU General Public License as published by |
| # the Free Software Foundation; either version 2 of the License, or |
| # (at your option) any later version. |
| # |
| # This program is distributed in the hope that it will be useful, |
| # but WITHOUT ANY WARRANTY; without even the implied warranty of |
| # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| # GNU General Public License for more details. |
| # |
| # A copy of the GNU General Public License is available at |
| # https://www.R-project.org/Licenses/ |
| |
| ARMAacf <- function(ar = numeric(), ma = numeric(), lag.max = r, |
| pacf = FALSE) |
| { |
| p <- length(ar) |
| q <- length(ma) |
| if(!p && !q) stop("empty model supplied") |
| r <- max(p, q + 1) |
| if(p > 0) { |
| if(r > 1) { |
| if(r > p) { ## pad with zeros so p >= q+1 |
| ar <- c(ar, rep(0, r - p)) |
| p <- r |
| } |
| p1 <- p + 1L |
| p2.1 <- p + p1 # = 2p + 1 |
| A <- matrix(0, p1, p2.1) |
| ind <- seq_len(p1) |
| ind <- as.matrix(expand.grid(ind, ind))[, 2L:1L] |
| ind[, 2] <- ind[, 1L] + ind[, 2L] - 1L |
| A[ind] <- c(1, -ar) |
| A[, 1L:p] <- A[, 1L:p] + A[, p2.1:(p + 2L)] |
| rhs <- c(1, rep(0, p)) |
| if(q > 0) { |
| psi <- c(1, ARMAtoMA(ar, ma, q)) |
| theta <- c(1, ma, rep(0, q+1L)) |
| for(k in 1L + 0:q) rhs[k] <- sum(psi * theta[k + 0:q]) |
| } |
| ind <- p1:1 |
| Acf <- solve(A[ind, ind], rhs) |
| Acf <- Acf[-1L]/Acf[1L] |
| } else Acf <- ar |
| if(lag.max > p) { |
| xx <- rep(0, lag.max - p) |
| Acf <- c(Acf, filter(xx, ar, "recursive", init = rev(Acf))) |
| } |
| Acf <- c(1, Acf[1L:lag.max]) |
| } else if(q > 0) { |
| x <- c(1, ma) |
| Acf <- filter(c(x, rep(0, q)), rev(x), sides=1)[-(1L:q)] |
| if(lag.max > q) Acf <- c(Acf, rep(0, lag.max - q)) |
| Acf <- Acf/Acf[1L] |
| } |
| names(Acf) <- 0:lag.max |
| if(pacf) drop(.Call(C_pacf1, Acf, lag.max)) else Acf |
| } |
| |
| acf2AR <- function(acf) |
| { |
| r <- as.double(drop(acf)) |
| order.max <- length(r) - 1 |
| if(order.max <= 0) stop("'acf' must be of length two or more") |
| z <- .Fortran(C_eureka, as.integer(order.max), r, r, |
| coefs = double(order.max^2), vars = double(order.max), |
| double(order.max)) |
| nm <- paste0("ar(",1L:order.max, ")") |
| matrix(z$coefs, order.max, order.max, dimnames=list(nm, 1L:order.max)) |
| } |
| |
| ARMAtoMA <- function(ar = numeric(), ma = numeric(), lag.max) |
| .Call(C_ARMAtoMA, as.double(ar), as.double(ma), as.integer(lag.max)) |