CODEWITHSUNDEEP
rarima
pdaawr
pdaawr

Reputation: 598

Save only one plot from tsdiag

I would like to retrieve and save to file only one plot out three from tsdiag function.

enter image description here

How can I achieve that? I've been trying to use additional parameter which but the plotting function is still returning the whole diagnosis.

tsdiag(arima_model, which = 1)  # does not work

I couldn't find anything in the documentation of stats either. It's quite easy to reproduce these plots manually but it would be great to get only one of them.

Upvotes: 1

Views: 642

Answers (2)

Georgi Boshnakov
Georgi Boshnakov

Reputation: 146

You could consider tsdiag.Sarima() from package sarima (https://cran.r-project.org/package=sarima), which provides the requested and further functionality, see https://geobosh.github.io/sarima/reference/tsdiag.Sarima.html. Another feature of tsdiag.Sarima() is that it uses the correct degrees of freedom to calculate p-values, which is another reason to use it. tsdiag.Sarima() offers also some alternative tests, see its help page or the link above.

For example, fit an ARIMA model to builtin AirPassengers data:

ap.arima <- arima(log(AirPassengers), order = c(0,1,1), seasonal = c(0,1,1))

This shows a graph similar to the default tsdiag one but with the d.f. as required by for the Ljung-Box test:

tsdiag.Sarima(ap.arima)

This shows also the Li-McLeod test (the residuals are dropped here):

tsdiag.Sarima(ap.arima, plot = 2:4)

Argument layout can be used to define a different layout of the plots, see ?layout for details. In the simplest case, it is a matrix with 1 column. This shows the autocorrelations and the LB p-values:

tsdiag.Sarima(ap.arima, plot = 2:3, layout = list(matrix(1:2, nrow = 2)))

plot can specify non-contiguous plots, of course. This shows the autocorrelations and the partial autocorrelations of the residuals

tsdiag.Sarima(ap.arima, plot = c(2,6), layout = list(matrix(1:2, nrow = 2)))

If argument plot specifies more plots than the layout can accommodate, you will get a menu to choose the desired plots:

tsdiag.Sarima(ap.arima, plot = 1:6, layout = list(matrix(1:2, nrow = 2)))

presents something like:

Select a plot number or 0 to exit 

1: residuals
2: acf of residuals
3: p values for Ljung-Box statistic
4: p values for Li-McLeod statistic
5: p values for Box-Pierce statistic
6: pacf of residuals

Selection:

layout allows you to allocate different amount of space to the plots, as well. This allocates 2/(1+2+2) of the vertical space (40%) to each of the second and third plot, and only 20% to the first (the standardised residuals):

tsdiag.Sarima(ap.arima, plot = c(1,2,6), 
    layout = list(matrix(1:3, nrow = 3),
    heights = c(1,2,2)))

The above works also for objects obtained from auto.arima. This modification of the example from the previous answer uses tsdiag.Sarima:

library("forecast")   # contains auto.arima
data(manaus, package = "boot")

m <- auto.arima(manaus)

tsdiag.Sarima(m)             # all plots
tsdiag.Sarima(m, plot = 2, layout = list(1))    # only 2nd
tsdiag.Sarima(m, plot = 2:3, layout = list(matrix(1:2, nrow = 2))) # plot 2 & 3

Upvotes: 1

tpetzoldt
tpetzoldt

Reputation: 5813

R is open source, so one can have a look in the sources and check if it is possible to create an own user-modified function.

The tsdiagfunction is found in package stats, source file ama0.R. Here a hacked version that supports a which-argument:

## modified from R package stats, file `ama0.R`
tsdiag.Arima <- tsdiag.arima0 <- function(object, gof.lag = 10, which = 1L:3L, ...) {
  ## plot standardized residuals, acf of residuals, Ljung-Box p-values
  oldpar <- par(mfrow = c(length(which), 1))
  on.exit(par(oldpar))
  rs <- object$residuals
  if (1L %in% which) {
    stdres <- rs/sqrt(object$sigma2)
    plot(stdres, type = "h", main = "Standardized Residuals", ylab = "")
    abline(h = 0)
  }
  if (2L %in% which) {
    acf(object$residuals, plot = TRUE, main = "ACF of Residuals",
        na.action = na.pass)
  }
  if (3L %in% which) {
    nlag <- gof.lag
    pval <- numeric(nlag)
    for(i in 1L:nlag) pval[i] <- Box.test(rs, i, type="Ljung-Box")$p.value
    plot(1L:nlag, pval, xlab = "lag", ylab = "p value", ylim = c(0,1),
         main = "p values for Ljung-Box statistic")
    abline(h = 0.05, lty = 2, col = "blue")
  }
}

Now we can test it:

library("forecast")   # contains auto.arima
library("boot")       # contains manaus data set

m <- auto.arima(manaus)

tsdiag(m)             # all plots
tsdiag(m, which=1)    # only 2nd
tsdiag(m, which=2:3)  # plot 2 and 3

Upvotes: 1

Related Questions