These tests are not actually computed in R. The seasonal package is
an interface to X-13-ARIMA-SEATS, which is where the computations are performed.
The QS statistics are based on a draft research report by Maravall (2012, Research Report, Banco de España). Details are also available in the documentation of X-13-ARIMA-SEATS.
These statistics are relatively straightforward to implement in R. Below I give you the general sketch to compute the statistic. (I give an example of the QS statistic computed for the original series).
# run the whole procedure for the "AirPassengers" series
require(seasonal)
m <- seas(AirPassengers)
The steps shown below follow the description given in the documentation of X-13-ARIMA-SEATS. First, the series for which the QS statistic is computed is differenced according to the chosen ARIMA model and the following rule:
$$
ndif = \max(1, \min(d + D, 2))
$$
where $ndif$ is the number of regular differences to be taken; $d$ and $D$ are respectively the number of regular and seasonal differences in the chosen ARIMA model. (If the QS statistic is computed for the series of residuals no differences are applied, $ndif=0$.)
require(polynom)
x <- AirPassengers
S <- frequency(x)
ndif <- max(1, min(sum(arimamodel(m)[c(2,5)]), 2))
dx <- filter(x, polynomial(c(1,-1))^ndif, sides=1)
dx <- window(dx, start=time(x)[ndif+1])
# alternatively, we can do it without package "polynom" by simply doing,
# if ndiff=1:
#dx <- diff(x)
# if ndiff=2:
#dx <- diff(diff(x))
Next, the first two autocorrelations of seasonal order (e.g. 12 and 24 in monthly data) are obtained. If these autocorrelations are lower or equal to zero, then they are set to zero.
R <- acf(dx, lag.max=S*2, plot=FALSE)$acf[-1,,1][c(S, 2*S)]
if (R[1] <= 0)
R[1] <- 0
if (R[2] <= 0)
R[2] <- 0
The statistic is defined as follows ($R_s$ and $R_{2s}$ denote the autocorrelations obtained in the previous step):
$$
QS = n(n+2)\left( \frac{R_s^2}{n-s} + \frac{R_{2s}^2}{n-2s} \right) \,,
$$
where $n$ is the number of observations in the differenced series and $s$ is
the periodicity of the data (12 in this case with monthly data). According to simulations exercises, this statistic follows approximately the $\chi^2$ distribuion with 2 degrees of freedom.
Thus, the statistic and the corresponding p-value can be obtained as follows:
n <- length(dx)
QS <- n*(n+2)*(R[1]^2/(n-S) + R[2]^2/(n-2*S))
pvalue <- pchisq(q=QS, df=2, lower.tail=FALSE)
round(c(QS=QS, p.value=pvalue), 4)
# QS p.value
# 167.6486 0.0000
which agrees with the output returned by seasonal::qs:
qs(m)["qsori",]
# qs p-val
# qsori 167.6486 0
In order to apply this test to other series, e.g. qsorievadj, you just need
to take this series from the output files and apply the operations shown above to it.
Given the large value of the QS statistic (and the implied low p-value), we can
conclude that there is seasonality in the series.
