# Computing the spectrum of a seasonal model in R

I am looking for a function which can:-

take a seasonal ARMA model

and return :-

the spectrum of the seasonal ARMA model.

I have seen:

library(astsa)
?arma.spec


However it's a bit clumsy to use arma.spec when my model has seasonal ar / ma terms.

For example (from the help page of arma.spec) we have that, we may compute the spectrum of a seasonal AR model like this:

Notice: How we have a seasonal AR term of order one, which is being fed to arma.spec by converting the seasonal lag to a non-seasonal lag:

arma.spec(ar=c(rep(0,11),.4), ma=.5, col=5, lwd=3, frequency=12)



I will illustrate my difficulty as follows:

> plot(AirPassengers)
> auto.arima(log(AirPassengers))
Series: log(AirPassengers)
ARIMA(0,1,1)(0,1,1)[12]

Coefficients:
ma1     sma1
-0.4018  -0.5569
s.e.   0.0896   0.0731

sigma^2 = 0.001371:  log likelihood = 244.7
AIC=-483.4   AICc=-483.21   BIC=-474.77


I do not know how to programmatically feed the above model to arma.spec.

I have tried:

auto.arima(log(AirPassengers))$coef arma.spec(auto.arima(log(AirPassengers))$coef,frequency=12)



I manually convert the seasonal MA to non-seasonal MA model and give it to arma.spec and I get a different picture.

arma.spec(ma=c(-.4,rep(0,10),-.55))



Hence I conclude that arma.spec does not understand seasonal models.

Query : How can I conveniently plot the spectrum of a seasonal ARMA process ?

Note : I wonder if there is a library which converts seasonal model to non-seasonal model which perhaps can then be used with arma.spec.

• It's been a long while since I've used ARIMA, but should the rep(0, 10) in your example be rep(0, 11) to align with the 13th month?
– Daniel V
Commented May 30 at 8:59

Here's a function that I think does what you want.

arma_model_spec <- function(object, ...) {
if (!inherits(object, "Arima")) {
stop("Argument must be an ARIMA model")
}
ar <- object$$model$$phi
ma <- object$$model$$theta
if (length(ar) == 0) ar <- 0
if (length(ma) == 0) ma <- 0
astsa::arma.spec(ar = ar, ma = ma, frequency = object$arma[5], ...) } forecast::auto.arima(log(AirPassengers)) |> arma_model_spec(main = "Spectrum of fitted model")  Created on 2024-06-05 with reprex v2.1.0 • Dear Sir, Many thanks for your reply. Commented Jun 9 at 6:36 Try this syntactic sugar on for size: extract_arma <- function(arima_model) { if (arima_model$arma[1] > 0) {
ar_input <- arima_model$coef[names(arima_model$coef)[grepl("^ar", names(arima_model$coef))]] } else { ar_input <- 0 } if (arima_model$arma[3] > 0) { ## Seasonal element of AR
ar_input <- c(ar_input,
rep(0, arima_model$arma[5] - arima_model$arma[1]),
arima_model$coef[names(arima_model$coef)[grepl("^sar", names(arima_model$coef))]]) } if (arima_model$arma[2] > 0) {
ma_input <- arima_model$coef[names(arima_model$coef)[grepl("^ma", names(arima_model$coef))]] } else { ma_input <- 0 } if (arima_model$arma[4] > 0) { ## Seasonal element of MA
ma_input <- c(ma_input,
rep(0, arima_model$arma[5] - arima_model$arma[2]),
arima_model$coef[names(arima_model$coef)[grepl("^sma", names(arima_model$coef))]]) } return(list(ar = ar_input, ma = ma_input freq = arima_model$arma[5]))
}


The function you've provided doesn't look like it's designed to receive a model as its inputs, only its values. As such, this function extracts the values in the format as close to what you've described as possible.

As I said in my comment, my memories of ARIMA are over a decade old so I can't remember exactly how to interface with it. My basic research suggests that the first seasonality term should start after the first complete cycle (i.e. after the 12th term, not on the 12th term).

The function can then be called as follows:

spec_input <- extract_arma(auto.arima(log(AirPassengers)))

arma.spec(ar = spec_input$ar, ma = spec_input$ma,
col = 5,
lwd = 3,
frequency = spec_input\$freq)


Hope this helps!

• Dear Daniel, Many thanks for your reply. I was wondering if there are R packages out there which can do this? To my knowledge there are two notations in place. 1. The Box-Jenkins and 2. the Signal processing notation for ARMA models. That is why it is safer to use a ready made package for this task, is what I feel. Commented Jun 3 at 8:38
• Not that I'm aware of - the packages that you're using I expect were developed independently with different design paradigms
– Daniel V
Commented Jun 4 at 4:15
• This answer, stats.stackexchange.com/questions/360133/…, talks about the notation in place for ARMA models in R. Commented Jun 9 at 5:21
• The forecast package has a function Arima and ?Arima says it is based on arima in base R. Hence I think the same notation is used by forecast package and base R. Commented Jun 9 at 5:56
• Also, the help of astsa:::sarima says in the Source that: This is an enhancement of ‘arima’ from the ‘stats’ package. Hence I think all 3 packages used in this conversation are using the same notation. Commented Jun 9 at 6:23