Distribution of residual autocorrelations for multiplicative seasonal ARMA models with uncorrelated but non-independent error terms

TL;DR

This paper extends portmanteau tests for SARMA models by relaxing the independence assumption on errors, analyzing residual autocorrelations under weaker conditions, and validating findings through simulations and real data application.

Contribution

It introduces new asymptotic results for residual autocorrelations in SARMA models with uncorrelated but dependent errors, broadening their applicability.

Findings

Asymptotic distributions of residual autocorrelations are derived.

Monte Carlo experiments validate the theoretical results.

Application to sunspot data demonstrates practical relevance.

Abstract

In this paper we consider portmanteau tests for testing the adequacy of multiplicative seasonal autoregressive moving-average (SARMA) models under the assumption that the errors are uncorrelated but not necessarily independent.We relax the standard independence assumption on the error term in order to extend the range of application of the SARMA models.We study the asymptotic distributions of residual and normalized residual empirical autocovariances and autocorrelations underweak assumptions on the noise. We establish the asymptotic behaviour of the proposed statistics. A set of Monte Carlo experiments and an application to monthly mean total sunspot number are presented.

Tables

Table 1. Table 1: Empirical size (in %) of the modified and standard versions of the LB and BP tests in the case of SARMA ( 0 , 1 ) ​ ( 0 , 1 ) s 0 1 subscript 0 1 𝑠 (0,1)(0,1)_{s} model ( 15 )-( 16 ) with α 1 = 0 . subscript 𝛼 1 0 \alpha_{1}=0. The nominal asymptotic level of the tests is α = 5 % 𝛼 percent 5 \alpha=5\% . The number of replications is N = 1 , 000 𝑁 1 000 N=1,000 .

$s$	Length $n$	Lag $m$	${\mathrm{LB}}_{\text{sn}}$	${\mathrm{BP}}_{\text{sn}}$	${\mathrm{LB}}_{\text{w}}$	${\mathrm{BP}}_{\text{w}}$	${\mathrm{LB}}_{\text{s}}$	${\mathrm{BP}}_{\text{s}}$
		$4$	4.3	4.0	6.8	6.8	8.7	8.7
		$8$	5.8	5.4	5.8	5.6	7.2	6.7
4	$n=500$	$12$	4.9	4.8	6.0	5.4	7.2	6.3
		$15$	6.3	5.3	5.1	4.7	5.8	5.3
		$18$	5.7	5.1	5.4	4.3	5.7	4.9
		$20$	5.7	5.0	5.2	3.8	5.4	4.1
		$4$	4.1	4.1	5.7	5.6	7.5	7.5
		$8$	4.2	4.2	4.7	4.7	5.9	5.6
4	$n=2,000$	$12$	4.9	4.8	3.9	3.9	4.6	4.4
		$15$	4.6	4.6	4.9	4.4	4.9	4.5
		$18$	4.0	3.8	4.4	4.2	4.6	4.1
		$20$	4.5	4.5	5.1	5.0	5.4	5.1
		$4$	4.7	4.5	12.8	12.1	20.8	20.6
		$8$	6.5	6.4	12.1	11.4	15.1	14.6
12	$n=500$	$12$	5.4	5.1	10.7	10.1	11.7	10.8
		$15$	6.0	5.6	8.7	8.4	9.8	9.4
		$18$	5.1	4.5	8.3	7.5	9.7	8.3
		$20$	5.5	4.3	8.6	7.6	9.8	8.4
		$4$	4.4	4.3	7.2	7.0	14.0	13.7
		$8$	4.9	4.9	5.8	5.7	9.9	9.8
12	$n=2,000$	$12$	4.5	4.4	6.0	5.6	7.0	6.6
		$15$	4.5	4.5	6.1	6.0	7.0	6.9
		$18$	5.4	5.3	6.0	5.7	6.9	6.5
		$20$	4.9	4.8	6.3	6.1	7.5	7.1

Table 2. Table 2: Empirical size (in %) of the modified and standard versions of the LB and BP tests in the case of SARMA ( 0 , 1 ) ​ ( 0 , 1 ) s 0 1 subscript 0 1 𝑠 (0,1)(0,1)_{s} model ( 15 )-( 16 ) with α 1 = 0.45 . subscript 𝛼 1 0.45 \alpha_{1}=0.45. The nominal asymptotic level of the tests is α = 5 % 𝛼 percent 5 \alpha=5\% . The number of replications is N = 1 , 000 𝑁 1 000 N=1,000 .

$s$	Length $n$	Lag $m$	${\mathrm{LB}}_{\text{sn}}$	${\mathrm{BP}}_{\text{sn}}$	${\mathrm{LB}}_{\text{w}}$	${\mathrm{BP}}_{\text{w}}$	${\mathrm{LB}}_{\text{s}}$	${\mathrm{BP}}_{\text{s}}$
		$4$	3.1	2.9	4.5	4.3	18.3	18.3
		$8$	2.8	2.6	3.3	3.2	11.7	11.3
4	$n=500$	$12$	3.5	3.2	2.8	2.4	9.7	9.3
		$15$	2.8	2.3	2.9	2.7	10.0	9.0
		$18$	1.9	1.6	2.7	2.4	9.1	7.8
		$20$	1.7	1.5	2.6	2.3	9.0	7.9
		$4$	5.5	5.4	4.2	4.2	19.9	19.9
		$8$	4.8	4.8	3.6	3.6	13.9	13.7
4	$n=2,000$	$12$	4.5	4.3	3.8	3.6	11.4	11.3
		$15$	4.8	4.8	3.6	3.6	10.7	10.5
		$18$	3.7	3.7	3.4	3.4	10.2	9.8
		$20$	3.2	3.1	3.5	3.5	9.4	9.1
		$4$	3.6	3.5	9.1	8.9	29.6	29.1
		$8$	4.0	3.8	6.4	6.0	20.4	19.7
12	$n=500$	$12$	4.2	3.5	5.5	5.2	15.4	14.6
		$15$	2.2	1.9	5.4	4.8	14.6	13.5
		$18$	1.7	1.4	4.3	4.2	13.3	12.0
		$20$	1.8	1.6	4.9	4.2	12.9	11.8
		$4$	4.4	4.4	4.5	4.5	27.1	27.1
		$8$	4.8	4.7	4.1	3.9	20.3	20.1
12	$n=2,000$	$12$	3.8	3.8	4.6	4.4	15.3	15.2
		$15$	4.6	4.5	3.5	3.4	13.9	13.7
		$18$	3.8	3.7	3.6	3.5	12.6	12.4
		$20$	3.5	3.3	4.1	3.8	11.9	11.8

Table 3. Table 3: Empirical power (in %) of the modified and standard versions of the LB and BP tests in the case of SARMA ( 0 , 1 ) ​ ( 0 , 1 ) s 0 1 subscript 0 1 𝑠 (0,1)(0,1)_{s} model ( 17 )-( 16 ) with α 1 = 0 subscript 𝛼 1 0 \alpha_{1}=0 . The number of replications is N = 1 , 000 𝑁 1 000 N=1,000 .

$s$	Length $n$	Lag $m$	${\mathrm{LB}}_{\text{sn}}$	${\mathrm{BP}}_{\text{sn}}$	${\mathrm{LB}}_{\text{w}}$	${\mathrm{BP}}_{\text{w}}$	${\mathrm{LB}}_{\text{s}}$	${\mathrm{BP}}_{\text{s}}$
		$4$	98.9	98.9	100.0	100.0	100.0	100.0
		$8$	98.1	98.0	100.0	100.0	100.0	100.0
4	$n=500$	$12$	96.3	96.3	100.0	100.0	100.0	100.0
		$15$	94.9	94.9	99.9	99.9	100.0	100.0
		$18$	93.7	93.5	99.8	99.7	100.0	100.0
		$20$	91.6	91.4	99.7	99.7	100.0	100.0
		$4$	100.0	100.0	100.0	100.0	100.0	100.0
		$8$	100.0	100.0	100.0	100.0	100.0	100.0
4	$n=2,000$	$12$	100.0	100.0	100.0	100.0	100.0	100.0
		$15$	100.0	100.0	100.0	100.0	100.0	100.0
		$18$	100.0	100.0	100.0	100.0	100.0	100.0
		$20$	100.0	100.0	100.0	100.0	100.0	100.0
		$4$	98.1	98.1	100.0	100.0	100.0	100.0
		$8$	97.7	97.7	99.9	99.9	100.0	100.0
12	$n=500$	$12$	96.6	96.6	99.9	99.9	100.0	100.0
		$15$	95.5	95.4	99.9	99.9	100.0	100.0
		$18$	92.9	92.7	99.9	99.9	100.0	100.0
		$20$	90.9	90.5	99.9	99.9	100.0	100.0
		$4$	100.0	100.0	100.0	100.0	100.0	100.0
		$8$	100.0	100.0	100.0	100.0	100.0	100.0
12	$n=2,000$	$12$	100.0	100.0	100.0	100.0	100.0	100.0
		$15$	100.0	100.0	100.0	100.0	100.0	100.0
		$18$	100.0	100.0	100.0	100.0	100.0	100.0
		$20$	100.0	100.0	100.0	100.0	100.0	100.0

Table 4. Table 4: Empirical power (in %) of the modified and standard versions of the LB and BP tests in the case of ARMA ( 0 , 1 ) ​ ( 0 , 1 ) s 0 1 subscript 0 1 𝑠 (0,1)(0,1)_{s} model ( 17 )-( 16 ) with α 1 = 0.45 subscript 𝛼 1 0.45 \alpha_{1}=0.45 . The number of replications is N = 1000 𝑁 1000 N=1000 .

$s$	Length $n$	Lag $m$	${\mathrm{LB}}_{\text{sn}}$	${\mathrm{BP}}_{\text{sn}}$	${\mathrm{LB}}_{\text{w}}$	${\mathrm{BP}}_{\text{w}}$	${\mathrm{LB}}_{\text{s}}$	${\mathrm{BP}}_{\text{s}}$
		$4$	89.2	89.1	98.5	98.5	100.0	100.0
		$8$	83.0	82.9	96.8	96.7	100.0	100.0
4	$n=500$	$12$	71.5	71.5	96.3	96.3	100.0	100.0
		$15$	64.9	64.7	95.6	95.6	100.0	100.0
		$18$	56.2	55.7	95.1	95.0	100.0	100.0
		$20$	50.6	49.4	94.5	94.4	100.0	100.0
		$4$	99.5	99.5	100.0	100.0	100.0	100.0
		$8$	99.4	99.4	100.0	100.0	100.0	100.0
4	$n=2,000$	$12$	99.1	99.1	100.0	100.0	100.0	100.0
		$15$	99.2	99.2	99.9	99.9	100.0	100.0
		$18$	99.5	99.5	100.0	100.0	100.0	100.0
		$20$	99.2	99.2	99.9	99.9	100.0	100.0
		$4$	90.1	90.1	98.4	98.3	100.0	100.0
		$8$	82.4	82.2	97.3	97.2	100.0	100.0
12	$n=500$	$12$	74.3	74.0	97.0	96.9	100.0	100.0
		$15$	68.7	68.3	96.9	96.5	100.0	100.0
		$18$	56.9	56.0	96.1	95.7	100.0	100.0
		$20$	48.9	47.9	96.0	95.7	100.0	100.0
		$4$	99.6	99.6	100.0	100.0	100.0	100.0
		$8$	99.6	99.6	99.9	99.9	100.0	100.0
12	$n=2,000$	$12$	99.5	99.5	100.0	100.0	100.0	100.0
		$15$	99.4	99.4	100.0	100.0	100.0	100.0
		$18$	99.5	99.5	100.0	100.0	100.0	100.0
		$20$	99.0	99.0	100.0	99.9	100.0	100.0

Equations

$$a_{{\theta }_0}(L){\mathbf{a}}_{{\theta }_0}(L)X_t=b_{{\theta }_0}(L){\mathbf{b}}_{{\theta }_0}(L){ϵ}_t,\forall t\in \mathbb{Z},$$

$$\left\{\text{strong SARMA}(p,q)(P,Q{)}_s\right\}\subset \left\{\text{weak SARMA}(p,q)(P,Q{)}_s\right\}.$$

$$Q_m^{\text{textsc}bp}=n\sum _{h=1}^m{\hat{\rho }}^2(h)\text{ and }{Q}_m^{\text{textsc}lb}=n(n+2)\sum _{h=1}^m\frac{{\hat{\rho }}^2(h)}{n-h},$$

$$\Theta$$

$$a_{\theta }(z)=1-\sum _{i=1}^p{a}_i{z}^i,{\mathbf{a}}_{\theta }(z)=1-\sum _{j=1}^P{\mathbf{a}}_j{z}^{sj},b_{\theta }(z)=1-\sum _{i=1}^p{b}_i{z}^i\text{ and }{\mathbf{b}}_{\theta }(z)=1-\sum _{j=1}^Q{\mathbf{b}}_j{z}^{sj}$$

$$\displaystyle\hskip 14.22636pt\text{ have all their zeros outside the unit disk and have no zero in common}\Big{\}}\ .$$

$${\alpha }_{ϵ}\left(h\right)=\underset{A\in {\mathcal{F}}_{-\infty }^t,B\in {\mathcal{F}}_{t+h}^{+\infty }}{\sup }\left|\mathbb{P}\left(A\cap B\right)-\mathbb{P}(A)\mathbb{P}(B)\right|,$$

$$e_t(\theta )$$

$$+\sum _{j=1}^Q{\mathbf{b}}_j{e}_{t-sj}(\theta )-\sum _{i=1}^q\sum _{j=1}^Q{b}_i{\mathbf{b}}_j{e}_{t-sj-i}(\theta ),$$

$${\mathrm{L}}_n(\theta ,{\sigma }^2)=\prod _{t=1}^n\frac{1}{(2\pi {)}^{1/2}\sqrt{{\sigma }^2}}\exp \left\{-\frac{e_t^2(\theta )}{2{\sigma }^2}\right\}.$$

$$({\hat{\theta }}_n,{\hat{\sigma }}^2)=\arg \underset{\theta ,{\sigma }^2}{\min }\left\{\log ({\sigma }^2)+\frac{1}{2{\sigma }^2}{Q}_n(\theta )\right\}\text{ where }{Q}_n(\theta )=\frac{1}{n}\sum _{t=1}^n{e}_t^2(\theta ).$$

$$\sqrt{n}\left({\hat{\theta }}_n-{\theta }_0\right)=-J^{-1}\frac{1}{\sqrt{n}}\sum _{t=1}^n{{\rm Y}}_t+{\mathrm{o}}_{\mathbb{P}}(1)\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}\mathcal{N}(0,\Sigma =J^{-1}I{J}^{-1}),$$

$$J=J({\theta }_0)=\frac{2}{{\sigma }_0^2}\mathbb{E}\left[\frac{\partial {ϵ}_t({\theta }_0)}{\partial \theta }\frac{\partial {ϵ}_t({\theta }_0)}{\partial {\theta }^{\prime }}\right],I=I({\theta }_0)=\sum _{h=-\infty }^{+\infty }\mathrm{Cov}({{\rm Y}}_t,{{\rm Y}}_{t-h})\text{ and }{{\rm Y}}_t=\frac{2}{{\sigma }_0^2}{ϵ}_t({\theta }_0)\frac{\partial {ϵ}_t({\theta }_0)}{\partial \theta }.$$

$${\hat{\gamma }}_m={\left(\hat{\gamma }(1),\dots ,\hat{\gamma }(m)\right)}^{\prime }\text{where}\hat{\gamma }(h)=\frac{1}{n}\sum _{t=h+1}^n{\hat{e}}_t{\hat{e}}_{t-h}\text{for}0\le h<n.$$

$${\hat{\rho }}_m={\left(\hat{\rho }(1),\dots ,\hat{\rho }(m)\right)}^{\prime }\text{where}\hat{\rho }(h)=\hat{\gamma }(h)/\hat{\gamma }(0).$$

$${\gamma }_m={\left(\gamma (1),\dots ,\gamma (m)\right)}^{\prime }\text{ where }\gamma (h)=\frac{1}{n}\sum _{t=h+1}^n{ϵ}_t{ϵ}_{t-h}\text{for}0\le h<n.$$

$$\displaystyle\Xi=\left(\begin{array}[]{cc}\Sigma&\Sigma_{\hat{\theta}_{n},\gamma_{m}}\\ \Sigma^{\prime}_{\hat{\theta}_{n},\gamma_{m}}&\Sigma_{\gamma_{m}}\end{array}\right)=\sum_{h=-\infty}^{+\infty}\mathbb{E}w_{t}w^{\prime}_{t-h},\text{ where }w_{t}=\left(\begin{array}[]{c}w_{1t}\\ w_{2t}\end{array}\right)\in\mathbb{R}^{k_{0}+m}$$

$$\sqrt{n}{\left({\hat{\theta }}_n^{\prime }-{\theta }_0^{\prime },{\gamma }_m^{\prime }\right)}^{\prime }=\frac{1}{\sqrt{n}}\sum _{t=1}^n{w}_t+{\mathrm{o}}_{\mathbb{P}}(1)\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}\mathcal{N}(0,\Xi ),$$

$$\hat{\gamma }(h)$$

$${\hat{\gamma }}_m:={\left(\hat{\gamma }(1),\dots ,\hat{\gamma }(m)\right)}^{\prime }={\gamma }_m+{\Phi }_m({\hat{\theta }}_n-{\theta }_0)+{\mathrm{O}}_{\mathbb{P}}(1/n),$$

$$\Phi_{m}=\mathbb{E}\left\{\left(\begin{array}[]{c}\epsilon_{t-1}\\ \vdots\\ \epsilon_{t-m}\end{array}\right)\frac{\partial\epsilon_{t}}{\partial\theta^{\prime}}\right\}.$$

$$\sqrt{n}{\hat{\gamma }}_m\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}\mathcal{N}\left(0,{\Sigma }_{{\hat{\gamma }}_m}\right)\text{mbox}and\sqrt{n}{\hat{\rho }}_m\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}\mathcal{N}\left(0,{\Sigma }_{{\hat{\rho }}_m}\right)\text{mbox}where,$$

$${\Sigma }_{{\hat{\gamma }}_m}$$

$${\hat{\Phi }}_m=\frac{1}{n}\sum _{t=1}^n\left\{{\left({\widehat{ϵ}}_{t-1},\dots ,{\widehat{ϵ}}_{t-m}\right)}^{\prime }\frac{\partial {\widehat{ϵ}}_t}{\partial {\theta }^{\prime }}\right\}\text{ and }{\hat{\sigma }}^2=\hat{\gamma }(0)=\frac{1}{n}\sum _{t=1}^n{\widehat{ϵ}}_t^2.$$

$$Z_m({\xi }_m)=\sum _{i=1}^m{\xi }_{i,m}{Z}_i^2$$

$$Q_m^{\text{textsc}lb}>S_m(1-\alpha )(\text{ resp. }{Q}_m^{\text{textsc}bp}>S_m(1-\alpha )),$$

$$\sqrt{n}{\hat{\gamma }}_m$$

$$C_m=\frac{1}{n^2}\sum _{t=1}^n{S}_t{S}_t^{\prime }\text{ where }{S}_t=\sum _{j=1}^t\left(\Lambda w_j-\Lambda \bar{w}\right)\text{ with }\bar{w}=\frac{1}{n}\sum _{t=1}^n{w}_t.$$

$${\mathcal{U}}_K=B_K^{\prime }(1)V_K^{-1}{B}_K(1)\text{ where }{V}_K={\int }_0^1\left(B_K(r)-r{B}_K(1)\right){\left(B_K(r)-r{B}_K(1)\right)}^{\prime }dr.$$

$$n{\hat{\gamma }}_m^{\prime }{C}_m^{-1}{\hat{\gamma }}_m\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}{\mathcal{U}}_m\text{ and }n{\sigma }_0^4{\hat{\rho }}_m^{\prime }{C}_m^{-1}{\hat{\rho }}_m\underset{n\to \infty }{\overset{\mathrm{d}}{\to }}{\mathcal{U}}_m$$

$$\hat{J}=\frac{2}{{\hat{\sigma }}^2}\frac{1}{n}\sum _{t=1}^n\frac{\partial e_t({\hat{\theta }}_n)}{\partial \theta }\frac{\partial e_t({\hat{\theta }}_n)}{\partial {\theta }^{\prime }}.$$

$$\hat{\Lambda}=\Big{(}\hat{\Phi}_{m}{|}I_{m}\Big{)}\text{ and }\hat{w}_{t}=\left(-2\hat{J}_{n}^{-1}\frac{\partial e_{t}(\hat{\theta}_{n})}{\partial\theta^{\prime}}\frac{1}{\hat{\sigma}^{2}}\hat{e}_{t},\hat{e}_{t}\hat{e}_{t-1},\dots,\hat{e}_{t}\hat{e}_{t-m}\right)^{\prime}.$$

Full text

62M10,

62F03,

62F05,

91B84,

62P05,

\startlocaldefs\endlocaldefs\doublespacing

and   t2Corresponding author

1 Introduction

The multiplicative seasonal autoregressive moving average (SARMA) model of order $(p,q)(P,Q)_{s}$ for the univariate time series $X=(X_{t})_{t\in\mathbb{Z}}$, is defined by

[TABLE]

where $\theta_{0}=(a_{01},\dots,a_{0p},b_{01},\dots,b_{0q},\mathbf{a}_{01},\dots,\mathbf{a}_{0P},\mathbf{b}_{01},\dots,\mathbf{b}_{0Q})^{\prime}$ and where the nonseasonal AR and MA operators are defined by $a_{\theta_{0}}(L)=1-\sum_{i=1}^{p}a_{0i}L^{i}$ and $b_{\theta_{0}}(L)=1-\sum_{i=1}^{q}b_{0i}L^{i}$, respectively, while the seasonal AR and MA operators are given by $\mathbf{a}_{\theta_{0}}(L)=1-\sum_{i=1}^{P}\mathbf{a}_{0i}L^{si}$ and $\mathbf{b}_{\theta_{0}}(L)=1-\sum_{i=1}^{Q}\mathbf{b}_{0i}L^{si}$, respectively, $s$ denotes the length of the seasonal period and $L$ stands for the backshift operator. It is assumed that the model defined by (1) is stationary, invertible and not redundant. Without loss of generality, we also assume that $a_{0p}^{2}+b_{0q}^{2}+\mathbf{a}_{0P}^{2}+\mathbf{b}_{0Q}^{2}\neq 0$ (by convention ${a}_{00}={b}_{00}=\mathbf{a}_{00}=\mathbf{b}_{00}=1$).

In the standard situation ${\epsilon}=({\epsilon}_{t})_{t\in\mathbb{Z}}$ is assumed to be a sequence of independent and identically distributed (iid for short) random variables with zero mean and common variance. In this standard framework, $(\epsilon_{t})$ is said to be a strong white noise and the representation (1) is called a strong SARMA$(p,q)(P,Q)_{s}$ process. In contrast with this previous definition, the representation (1) is said to be a weak SARMA$(p,q)(P,Q)_{s}$ if the noise process $(\epsilon_{t})$ is a weak white noise, that is, if it satisfies

• (A0):

$\mathbb{E}({{\epsilon}}_{t})=0$, $\mbox{Var}\left({{\epsilon}}_{t}\right)=\sigma^{2}_{0}$ and $\mbox{Cov}\left({{\epsilon}}_{t},{{\epsilon}}_{t-h}\right)=0$ for all $t\in\mathbb{Z}$ and all $h\neq 0$.

A strong white noise is obviously a weak white noise, because independence entails uncorrelatedness, but the reverse is not true. It is clear from these definitions that the following inclusions hold:

[TABLE]

After estimating the SARMA process, the next important step in the modeling consists in checking if the estimated model fits satisfactorily the data. Thus, under the null hypothesis that the model has been correctly identified, the residuals ($\hat{\epsilon}_{t}$) are approximately a white noise. This adequacy checking step validates or invalidates the choice of the orders $(p,q)$ and $(P,Q)_{s}$.

Based on the residual empirical autocorrelations $\hat{\rho}(h)=\sum_{t=1+h}^{n}\hat{\epsilon}_{t}\hat{\epsilon}_{t-h}/\sum_{t=1}^{n}\hat{\epsilon}_{t}^{2}$, where $n$ is the length of the series, [8] have proposed a goodness-of-fit test, the so-called portmanteau test, for strong ARMA models. A modification of their test has been proposed by [21]. It is nowadays one of the most popular diagnostic checking tools in ARMA modeling of time series. These tests are defined by

[TABLE]

where $m$ is a fixed integer. The statistic ${Q}_{m}^{\textsc{lb}}$ has the same asymptotic chi-squared distribution as $Q_{m}^{\mathrm{BP}}$ and has the reputation of doing better for small or medium sized sample (see [21]). For weak ARMA models, [13] show that the asymptotic distributions of the statistics defined in (2) are no longer chi-square distributions but a mixture of chi-squared distributions, weighted by eigenvalues of the asymptotic covariance matrix of the vector of autocorrelations. Recently, [7] proposed an alternative method based on a self-normalization approach to construct a new test statistic which is asymptotically distribution-free under the null hypothesis.

In many situations, these tests are implemented to check the lack of fit of SARMA models. However, the traditional methodology of Box and Jenkins cannot be extended to the case of SARMA models when $s>1$ because of the multiplicative contraints on the parameters. This standard methodology needs to be adapted to take into account the possible lack of independence of the errors terms. See, for instance, [12] and [23] who considered serial correlation testing in multiplicative seasonal univariate time series models. Duchesne (see [12]) proposed his test statistic based on a kernel-based spectral density estimator, whose weighting scheme is more adapted to autocorrelations associated to seasonal lags. The standard tests procedure consist in rejecting the null hypothesis of a SARMA$(p,q)(P,Q)_{s}$ model if the statistics (2) are larger than a certain quantile of a chi-squared distribution with $m-(p+q+P+Q)>0$ degrees of freedom. Consequently, these standard tests are not applicable for $m\leq p+q+P+Q$.

The works on the portmanteau statistic of SARMA models are generally performed under the assumption that the errors $\epsilon_{t}$ are independent. This independence assumption is often considered too restrictive by practitioners. It precludes conditional heteroscedasticity and/or other forms of nonlinearity (see [17], for a review on weak univariate ARMA models). In this framework, we relax the standard independence assumption on the error term in order to be able to cover SARMA representations of general nonlinear models. For the asymptotic theory of weak SARMA, notable exception are [6] where the consistency and the asymptotic normality of the quasi-maximum likelihood estimator (QMLE) for weak multivariate SARMA models are studied. They also study a particular case of the asymptotic distributions of residual autocovariances and autocorrelations at the seasonal lags $1s,2s,3s,\dots,ms$.

This paper is devoted to the problem of the validation step of weak SARMA representations. We consider portmanteau test statistics based on the residual empirical autocorrelations but not necessarily at multiple lags $s$ as in [6]. For such models, we show that the asymptotic distributions of the statistics defined in (2) are no longer chi-square distributions but a mixture of chi-squared distributions, weighted by eigenvalues of the asymptotic covariance matrix of the vector of autocorrelations as in [6]. We also proposed another modified statistics based on a self-normalization approach which are asymptotically distribution-free under the null hypothesis and generalize the result of [7].

In Monte Carlo experiments, we illustrate that the proposed test statistics have reasonable finite sample performance. Under nonindependent errors, it appears that the standard test statistics are generally non reliable, overrejecting severely, while the proposed tests statistics offer satisfactory levels. Even for independent errors, they seem preferable to the standard ones, when the number $m$ of autocorrelations is small. Moreover, the error of first kind is well controlled. Contrarily to the standard tests (2), the proposed tests can be used safely for $m$ small (see for instance Figure 1). For all these above reasons, we think that the modified versions that we propose in this paper are preferable to the standard ones for diagnosing SARMA models under nonindependent errors. Other contribution is to improve the results concerning the statistical analysis of weak SARMA models by considering the adequacy problem.

The article is organized as follows. In the next section, we briefly recall the results on the QMLE asymptotic distribution obtained by [6] when $(\epsilon_{t})$ satisfies mild mixing assumptions. We study the asymptotic behaviour of the residuals autocovariances and autocorrelations under weak assumptions on the noise in Section 3.1. It is also shown how the standard portmanteau tests (2) must be adapted in the case of multiplicative seasonal ARMA models with non-independent innovations. In Section 3.2 we derive the asymptotic distribution of residuals autocovariances and autocorrelations using self-normalization approach and we establish the asymptotic behaviour of the proposed statistics. Section 4 proposes numerical illustrations and an illustrative application on real data. We provide a conclusion in Section 5. The technical proofs are relegated to the appendix.

2 Estimating weak SARMA models

In this section, we recall the results on the QMLE asymptotic distribution obtained by [6] when $(\epsilon_{t})$ satisfies mild mixing assumptions in order to have a self-containing paper.

The unknown parameter of interest $\theta_{0}$ is supposed to belong to the parameter space

[TABLE]

To ensure the asymptotic theory of the QMLE, we assume that the parametrization satisfies the following smoothness conditions. Without loss of generality, we may assume that $\Theta$ is compact.

• (A1):

The process $\epsilon=(\epsilon_{t})_{t\in\mathbb{Z}}$ is ergodic and strictly stationary.

For the asymptotic normality of the QMLE, additional assumptions are required. It is necessary to assume that $\theta_{0}$ is not on the boundary of the parameter space ${\Theta}$.

• (A2):

We have $\theta_{0}\in\stackrel{{\scriptstyle\circ}}{{{\Theta}}}$, where $\stackrel{{\scriptstyle\circ}}{{{\Theta}}}$ denotes the interior of ${\Theta}$.

To control the serial dependence of the stationary process $(\epsilon_{t})$, we introduce the strong mixing coefficients $\alpha_{\epsilon}(h)$ defined by

[TABLE]

where $\mathcal{F}_{-\infty}^{t}=\sigma(\epsilon_{u},u\leq t)$ and $\mathcal{F}_{t+h}^{+\infty}=\sigma(\epsilon_{u},u\geq t+h)$. We use $|\cdot|$ to denote the Euclidian norm $|z|=\left(\sum_{i=1}^{\tau}z_{i}^{2}\right)^{1/2}$ of a column vector $z=(z_{1},\dots,z_{\tau})^{\prime}$. We will make an integrability assumption on the moment of the noise and a summability condition on the strong mixing coefficients $(\alpha_{{\epsilon}}(k))_{k\geq 0}$.

• (A3):

\text{We have }\mathbb{E}\big{|}{\epsilon}_{t}|^{4+2\nu}<\infty\text{ and }\sum_{k=0}^{\infty}\left\{\alpha_{{\epsilon}}(k)\right\}^{\frac{\nu}{2+\nu}}<\infty\text{ for some }\nu>0.

Assumption (A3) from [14, 17] is a technical condition for proving the asymptotic theory of the QMLE. The integrability assumption on the moment of the noise is not very restrictive in this framework because the innovation process $(\epsilon_{t})$ is directly observed (see [25]).

For the estimation of SARMA and multivariate SARMA models, the commonly used estimation method is the quasi-maximum likelihood estimation, which can be also viewed as a nonlinear least squares estimation (LSE). Given a realization ${X}_{1},{X}_{2},\dots,{X}_{n}$ satisfying (1), the variable $\epsilon_{t}(\theta)$ can be approximated, for $0<t\leq n,$ by ${e}_{t}(\theta)$ defined recursively by

[TABLE]

where the unknown initial values are set to zero: ${e}_{0}(\theta)=\dots={e}_{1-q-sQ}(\theta)={X}_{0}=\dots={X}_{1-p-sP}=0$. The Gaussian quasi-likelihood is given by

[TABLE]

A QMLE of $(\theta,\sigma^{2})$ is a measurable solution $(\hat{\theta}_{n},\hat{\sigma}^{2})$ of

[TABLE]

In all the sequel, we denote by $\xrightarrow[]{\mathrm{d}}$, the convergence in distribution. The symbol $\mathrm{o}_{\mathbb{P}}(1)$ is used for a sequence of random variables that converges to zero in probability. Under the above assumptions, [6] showed that $\hat{\theta}_{n}\to\theta_{0}\;a.s.$ as $n\to\infty$ and

[TABLE]

where

[TABLE]

Note that, the existence of the matrix $I(\theta_{0})$ is a consequence of (A3) and of Davydov’s inequality [10].

3 Diagnostic checking in weak SARMA models

In order to check the validity of the SARMA$(p,q)(P,Q)_{s}$ model, it is a common practice to examine the QMLE residuals $\hat{\epsilon}_{t}=\hat{e}_{t}={e}_{t}(\hat{\theta}_{n})$ where ${e}_{t}(\theta)$ is given by (2) for all $\theta\in\mathbb{R}^{k_{0}}$. For a fixed integer $m\geq 1$, consider the vector of residual autocovariances

[TABLE]

In the sequel, we will also need the vector of the first $m$ sample autocorrelations

[TABLE]

The statistics (2) are usually used to test the following null hypothesis

• ($\mathbf{H0}$) :

$({X}_{t})_{t\in\mathbb{Z}}$ satisfies a SARMA$(p,q)(P,Q)_{s}$ representation;

against the alternative

• ($\mathbf{H1}$) :

$({X}_{t})_{t\in\mathbb{Z}}$ does not admit a SARMA$(p,q)(P,Q)_{s}$ representation or $({X}_{t})_{t\in\mathbb{Z}}$ satisfies a SARMA$(p^{\prime},q^{\prime})(P^{\prime},Q^{\prime})_{s}$ representation with $p^{\prime}>p$ or $q^{\prime}>q$ or $P^{\prime}>P$ or $Q^{\prime}>Q$.

3.1 Asymptotic distribution of the residual autocorrelations

First note that the mixing assumptions (A3) entail the asymptotic normality of the "empirical" autocovariances

[TABLE]

It should be noted that $\gamma(h)$ is not a computable statistic because it depends on the unobserved innovations $\epsilon_{t}=\epsilon_{t}(\theta_{0})$ except when $p=q=P=Q=0$. Define the matrix

[TABLE]

with $w_{1t}=-J^{-1}\Upsilon_{t}=-2\sigma_{0}^{-2}\epsilon_{t}J^{-1}\left(\partial\epsilon_{t}/\partial\theta\right)$ and $w_{2t}=\left(\epsilon_{t-1},\dots,\epsilon_{t-m}\right)^{\prime}\epsilon_{t}$. Note that, the existence of the matrix $\Xi$ is a consequence of (A3) and of Davydov’s inequality [10].

The asymptotic distribution of $\sqrt{n}\hat{\rho}_{m}$ will be obtained from the joint asymptotic distribution of

[TABLE]

by applying the central limit theorem for mixing processes (see [19]).

Now, considering $\hat{\gamma}(h)$ and $\gamma(h)$ as values of the same function at the points $\hat{\theta}_{n}$ and $\theta_{0}$, a Taylor expansion about $\theta_{0}$ gives

[TABLE]

where $\theta_{n}^{*}$ is between $\hat{\theta}_{n}$ and $\theta_{0}.$ The last equality follows from the consistency of $\hat{\theta}_{n}$ and the fact that $\left(\partial\epsilon_{t-h}/\partial\theta^{\prime}\right)(\theta_{0})$ is not correlated with $\epsilon_{t}$ when $h\geq 0.$ Then for $h=1,\dots,m,$ we have

[TABLE]

where

[TABLE]

The following Proposition, which is a generalization of Theorem 4 obtained by [6], gives the limiting distribution of the residual autocovariances and autocorrelations of SARMA models.

Proposition 1**.**

When, $p>0$, $P>0$, $q>0$ and $Q>0$, under the above assumptions, we have

[TABLE]

[TABLE]

The proof of this result is similar to that given by [6] for Theorem 4.

The asymptotic variance matrices $\Sigma_{\hat{\gamma}_{m}}$ and $\Sigma_{\hat{\rho}_{m}}$ depend on the unknown matrices $\Xi$, $\Phi_{m}$ and the scalar $\sigma_{0}^{2}$. Matrix $\Phi_{m}$ and $\sigma_{0}^{2}$ can be estimated by its empirical counterpart, respectively

[TABLE]

Note that the matrix $(2\pi)^{-1}\Xi$ is the spectral density at frequency zero of the process $(w_{t})$, thus an estimator of $\Xi$ is given in Theorem 6 of [6]. Other estimators of such long-run variances are available in the literature (see for instance [1], [2], [11], [24], for general references). For the numerical illustrations presented in this paper, we used a Vector AR (VAR) spectral estimator given in Theorem 6 of [6] consisting in: i) fitting VAR$(r)$ models for $r=1,\dots,r_{\max}$ to the series $\hat{w}_{t}$, $t=1,\dots,n$, where $\hat{w}_{t}$ is obtained by replacing $\theta_{0}$ by $\hat{\theta}_{n}$ in $w_{t}$; ii) selecting the order $r$ which minimizes an information criterion and approximating $\Xi$ by $(2\pi)$ times the spectral density at frequency zero of the estimated VAR$(r)$ model. Hereafter, we used the AIC model selection criterion with $r_{\max}=5$.

From Proposition 1 we can deduce the following result, which gives the exact limiting distribution of the standard portmanteau statistics (2) under general assumptions on the innovation process of the fitted SARMA$(p,q)(P,Q)_{s}$ model.

Theorem 2**.**

Under Assumptions in Proposition 1 and ${\bf(H0)}$, the statistics $Q_{m}^{\textsc{lb}}$ and $Q_{m}^{\textsc{bp}}$ converge in distribution, as $n\rightarrow\infty,$ to

[TABLE]

where $\xi_{m}=(\xi_{1,m},\dots,\xi_{m,m})^{\prime}$ is the vector of the eigenvalues of the matrix $\Sigma_{\hat{\rho}_{m}}=\sigma_{0}^{-4}\Sigma_{\hat{\gamma}_{m}}$ and $Z_{1},\dots,Z_{m}$ are independent $\mathcal{N}(0,1)$ variables.

As in [6], Theorem 2 shows that for the asymptotic distribution of $Q_{m}^{\textsc{lb}}$ and $Q_{m}^{\textsc{bp}}$, the $\chi_{m-k_{0}}^{2}$ approximation is no longer valid in the framework of weak SARMA$(p,q)(P,Q)_{s}$ models. The true asymptotic distribution depends on nuisance parameters involving $\sigma_{0}^{2}$, the matrix $\Phi_{m}$ and the elements of $\Xi$. Consequently, in order to obtain the asymptotic distribution of the portmanteau statistics (2) under weak assumptions on the noise, one needs a consistent estimator of the asymptotic covariance matrix $\Sigma_{\hat{\rho}_{m}}$. We let $\hat{\Sigma}_{\hat{\rho}_{m}}$ the matrix obtained by replacing $\Xi$ by $\hat{\Xi}$, $\Phi_{m}$ by $\hat{\Phi}_{m}$ and $\sigma_{0}^{2}$ by $\hat{\sigma}^{2}$ in $\Sigma_{\hat{\rho}_{m}}$. Denote by $\hat{\xi}_{m}=(\hat{\xi}_{1,m},\dots,\hat{\xi}_{m,m})^{\prime}$ the vector of the eigenvalues of $\hat{\Sigma}_{\hat{\rho}_{m}}$. At the asymptotic level $\alpha$, the LB (Ljung-Box) test (resp. the BP (Box-Pierce) test) consists in rejecting the adequacy of the weak SARMA$(p,q)(P,Q)_{s}$ model when

[TABLE]

where $S_{m}(1-\alpha)$ is such that $\mathbb{P}\left\{Z_{m}(\hat{\xi}_{m})>S_{m}(1-\alpha)\right\}=\alpha$. We emphasize the fact that the proposed modified versions of the Box-Pierce and Ljung-Box statistics are more difficult to implement because their critical values have to be computed from the data.

3.2 Self-normalized asymptotic distribution of the residual autocorrelations

The nonparametric kernel estimator (see [1, 24]), used to estimate the matrix $\Xi$ causes serious difficulties regarding the choice of the sequence of weights. The parametric approach based on an autoregressive estimate of the spectral density of $w_{t}$ studied for instance by [2, 3, 4, 6, 11] is also facing the problem of choosing the truncation parameter. So the choice of the order of truncation is often crucial and difficult. In this section, we propose as in [7] an alternative method where we do not estimate an asymptotic covariance matrix. It is based on a self-normalization based approach to construct a test-statistic which is asymptotically distribution-free under the null hypothesis (see [7], for a reference in the ARMA cases). The idea comes from [22] and has been already extended by [20, 26, 27, 28] to more general frameworks. See also [29] for a review on some recent developments on the inference of time series data using the self-normalized approach. In this case, the critical values are not computed from the data since they are tabulated. In some sense, this method is finally closer to the standard method in which the critical values are simply deduced from a $\chi^{2}$-table.

We denote $\Lambda$ the matrix in $\mathbb{R}^{m\times(k_{0}+m)}$ defined in block formed by ${\Lambda}=({\Phi}_{m}{|}I_{m})$, where $I_{m}$ is the identity matrix of order $m$. In view of (4) and (10), we deduce that

[TABLE]

Contrarily to Subsection 3.1, we do not rely on the classical method that would consist in estimating the asymptotic covariance matrix of $\Lambda w_{t}$. We need to apply the functional central limit theorem holds for the process $w=(w_{t})_{t\geq 1}$ (see Lemma 1 in [22]).

Finally, we define the normalization matrix $C_{m}\in\mathbb{R}^{m\times m}$ by

[TABLE]

To ensure the invertibility of the normalization matrix $C_{m}$ which is proved in Lemma 6 of [7], we need the following technical assumption on the distribution of $\epsilon_{t}$.

• (A4):

The process $(\epsilon_{t})_{t\in\mathbb{Z}}$ has a positive density on some neighborhood of zero.

Let $(B_{K}(r))_{r\geq 0}$ be a $K$-dimensional Brownian motion starting from [math]. For $K\geq 1$, we denote $\mathcal{U}_{K}$ the random variable defined by

[TABLE]

The following theorem states the asymptotic distributions of the sample autocovariances and autocorrelations.

Theorem 3**.**

We assume that $p>0$, $q>0$, $P>0$ or $Q>0$. Under Assumptions of Proposition 1, (A4) and under the null hypothesis (H0), we have

[TABLE]

The proof of this result is postponed to Section A.

Of course, the above theorem is useless for practical purpose, because it does not involve any observable quantities. In practice, one has to replace the matrix $C_{m}$ and the variance of the noise $\sigma_{0}^{2}$ by their empirical or observable counterparts. The matrix $J$ can be easily estimated by his empirical counterpart

[TABLE]

Thus we define

[TABLE]

Finally we denote the normalization matrix $\hat{C}_{m}\in\mathbb{R}^{m\times m}$ by

[TABLE]

The following result is the applicable counterpart of Theorem 3.

Theorem 4**.**

Assume that $p>0$, $q>0$, $P>0$ or $Q>0$. Under Assumptions of Theorem 3, we have

[TABLE]

The proof of this result is postponed to Section A.

Based on the above result, we propose a modified version of the Ljung-Box statistic when one uses the statistic

[TABLE]

where the matrix $D_{n,m}\in\mathbb{R}^{m\times m}$ is diagonal with $((n+2)/(n-1),...,(n+2)/(n-m))$ as diagonal terms.

4 Numerical illustrations

In this section, by means of Monte Carlo experiments, we investigate the finite sample properties of the modified version of the portmanteau tests that we introduced in this work. The numerical illustrations of this section are made with the open source statistical software R (see R Development Core Team, 2017) or (see http://cran.r-project.org/).

4.1 Simulated models

First of all, we introduce the models that we simulate and we indicate the conventions that we adopt in the discussion and in the tables:

• •

$\mathrm{{LB}_{\textsc{w}}}$ and $\mathrm{BP_{\textsc{w}}}$ refer to modified LB and BP tests using ${Q}_{m}^{\textsc{lb}}$ and $Q_{m}^{\textsc{bp}}$ in Section 3.1

• •

$\mathrm{{LB}_{\textsc{s}}}$ and $\mathrm{BP_{\textsc{s}}}$ refer to LB and BP tests using the standard statistics (2).

• •

$\mathrm{{{LB}}_{\textsc{sn}}}$ and $\mathrm{{{BP}}_{\textsc{sn}}}$ refer to modified tests using the self-normalized statistics in Section 3.2

To generate the strong and the weak SARMA models, we consider the following SARMA$(0,1)(0,1)_{s}$ model

[TABLE]

with $\theta_{0}=(b_{01},\mathbf{b}_{01})^{\prime}=(-0.6,-0.7)^{\prime}$ and the innovation process $(\epsilon_{t})$ follows a strong or weak white noise.

The generalized autoregressive conditional heteroscedastic (GARCH) models is an important example of weak white noises in the univariate case (see [18]). So we first assume that in (15) the innovation process $\epsilon$ is the following ARCH$(1)$ model defined by

[TABLE]

where $(\eta_{t})_{t\geq 1}$ is a sequence of iid standard Gaussian random variables. To generate the strong SARMA, we assume that in (15) the innovation process follows (16) with $\alpha_{1}=0$.

4.2 Empirical size

We first simulate $N=1,000$ independent trajectories of size $n=2,000$ of models (15). The same series is partitioned as two series of sizes $n=500$ and $n=2,000$. For each of these $N$ replications, we use the quasi-maximum likelihood estimation method to estimate the coefficient $\theta_{0}$ and we apply portmanteau tests to the residuals for different values of $m\in\{4,8,12,15,18,20\}$, where $m$ is the number of autocorrelations used in the portmanteau test statistic. For the nominal level $\alpha=5\%$, the empirical size over the $N$ independent replications should vary between the significant limits 3.6% and 6.4% with probability 95% and belong to $[3.2\%,6.9\%]$ with a probability 99%. When the relative rejection frequencies are outside the 95% significant limits, they are displayed in bold type and they are underlined when they are outside the 99% significant limits in Tables 1 and 2.

For the standard Box-Pierce test, the model is therefore rejected when the statistic $Q_{m}^{\textsc{bp}}$ or $Q_{m}^{\textsc{lb}}$ is larger than $\chi_{(m-2)}^{2}(0.95)$ in a SARMA$(0,1)(0,1)_{s}$ case (see [23]). Consequently the empirical size is not available (n.a.) for the statistic $Q_{m}^{\textsc{bp}}$ or $Q_{m}^{\textsc{lb}}$ because they are not applicable for $m\leq 2$. For the proposed self-normalized test $\mathrm{BP}_{\textsc{sn}}$ or $\mathrm{LB}_{\textsc{sn}}$, the model is rejected when the statistic $Q_{m}^{\textsc{sn}}$ or $\tilde{Q}_{m}^{\textsc{sn}}$ is larger than $\mathcal{U}_{m}(0.95)$, where the critical values $\mathcal{U}_{K}(0.95)$ (for $K=1,\dots,20$) are tabulated in Lobato (see Table 1 in [22]).

Table 1 displays the relative rejection frequencies of the null hypothesis ${\bf(H0)}$ that the data generating process (DGP for short) follows a strong SARMA model (15)–(16) with $\alpha_{1}=0$, over the $N$ independent replications. When the seasonal period is $s=4$, for all tests, the percentages of rejection belong to the confident interval with probabilities 95% and 99%, except for $\mathrm{{LB}}_{\textsc{s}}$ and ${\mathrm{BP}}_{\textsc{s}}$ when $m=4$. Consequently all these tests well control the error of first kind. In contrast, when $s=12$, our proposed tests well also control the error of first kind (except for $\mathrm{{LB}}_{\textsc{w}}$ and ${\mathrm{BP}}_{\textsc{w}}$ when $n=500$) contrarily to the standard tests $\mathrm{{LB}}_{\textsc{s}}$, ${\mathrm{BP}}_{\textsc{s}}$ for all sizes. We draw the conclusion that, in this strong SARMA case, the proposed modified version may be clearly preferable to the standard ones.

Now, we repeat the same experiments on a weak SARMA models. As expected, Table 2 shows that the standard $\mathrm{{LB}}_{\textsc{s}}$ or $\mathrm{{BP}}_{\textsc{s}}$ test poorly performs in assessing the adequacy of this particular weak SARMA model. It can be seen that: 1) the observed relative rejection frequencies of $\mathrm{{LB}}_{\textsc{s}}$ and $\mathrm{{BP}}_{\textsc{s}}$ are definitely outside the significant limits, 2) the errors of the first kind are only globally well controlled by the proposed tests, for all $s$ when $n$ is large. We also tried the case where the ARCH$(1)$ model (16) have infinite fourth moments. As showing in Figure 1, the results are qualitatively similar to what we observe here.

Figure 1 displays the residual autocorrelations of a realization of size $n=2,000$ for weak SARMA$(0,1)(0,1)_{12}$ model (15)-(16) with $\alpha_{1}=1.3$ and their 5% significance limits under the strong SARMA$(0,1)(0,1)_{12}$ and weak SARMA$(0,1)(0,1)_{12}$ assumptions. This figure confirms clearly the conclusions drawn from Table 2. The horizontal dotted lines (blue color) correspond to the 5% significant limits obtained under the strong SARMA assumption. The solid lines (red color) and dashed lines (green color) correspond also to the 5% significant limits under the weak SARMA assumption. The full lines correspond to the asymptotic significance limits for the residual autocorrelations obtained in Proposition 1. The dashed lines (green color) correspond to the self-normalized asymptotic significance limits for the residual autocorrelations obtained in Theorem 4.

In these Monte Carlo experiments, we illustrate that the proposed test statistics have reasonable finite sample performance. Under nonindependent errors, it appears that the standard test statistics are generally non reliable, overrejecting severely, while the proposed tests statistics offer satisfactory levels. Even for independent errors, they seem preferable to the standard ones, when the number $m$ of autocorrelations is small and when $s=12$. Moreover, the error of first kind is well controlled. Contrarily to the standard tests based on $\mathrm{{BP}}_{\textsc{s}}$ or $\mathrm{{{LB}}_{\textsc{s}}}$, the proposed tests can be used safely for $m$ small (see for instance Figure 1). For all these above reasons, we think that the modified versions that we propose in this paper are preferable to the standard ones for diagnosing SARMA models under nonindependent errors.

4.3 Empirical power

In this section we repeat the same experiments as in Section 4.1 to examine the power of the tests for the null hypothesis of a SARMA$(0,1)(0,1)_{s}$ against the following SARMA$(1,1)(0,1)_{s}$ alternative defined by

[TABLE]

with $\theta_{0}=(a_{01},b_{01},\mathbf{b}_{01})^{\prime}=(0.8,-0.6,-0.7)^{\prime}$ and where the innovation process $\epsilon$ follows a strong or weak white noise introduced in Section 4.1. For each of these $N$ replications we fit a SARMA$(0,1)(0,1)_{s}$ models and perform standard and modified tests based on $m=4,8,12,15$, $18$ and $20$ residual autocorrelations.

Tables 3 and 4 compare the empirical powers of Model (17)-(16) with $\alpha_{1}=0$ and $\alpha_{1}=0.45$ respectively over the $N$ independent replications. For these particular strong and weak SARMA models, we notice that the standard $\mathrm{{BP}}_{\textsc{s}}$ and $\mathrm{{{LB}}_{\textsc{s}}}$ and our proposed tests have very similar powers except for $\mathrm{{BP}}_{\textsc{sn}}$ and $\mathrm{{{LB}}_{\textsc{sn}}}$ when $n=500$ in the weak case.

4.4 Application to real data

We now consider an application to monthly mean total sunspot number obtained by taking a simple arithmetic mean of the daily total sunspot number over all days of each calendar month. The observations (sunspot) covered the period from January 01, 2010 to December 31, 2018 which correspond to $n=108$ observations. The series exhibit seasonal behavior $(s=12)$. The data were obtain from the website of the World Data Center, Solar Influences Data Analysis Center, Royal Observatory of Belgium (http://www.sidc.be/silso/datafiles).

Let $Z_{t}=\log(\mbox{sunspot}_{t})-\log(\mbox{sunspot}_{t-1})$ and denoting by $X_{t}=Z_{t}-\mathbb{E}(Z_{t})$ the mean-corrected series. We adjust the particular SARMA$(3,1)(0,1)_{12}$ model of the form

[TABLE]

The quasi-maximum likelihood estimators of $\theta_{0}=(a_{03},\mathbf{b}_{01},\mathbf{b}_{01})^{\prime}$ were obtained as

[TABLE]

where the estimated asymptotic standard errors obtained from (4) (respectively the $p$-values), of the estimated parameters (first column), are given into brackets (respectively in parentheses). We apply portmanteau tests to the residuals of this model. Figure 2 displays the residual autocorrelations and their 5% significance limits under the strong SARMA$(3,1)(0,1)_{12}$ and weak SARMA$(3,1)(0,1)_{12}$ assumptions. In view of Figure 2, the diagnostic checking of residuals does not indicate any inadequacy. All of the sample autocorrelations should lie between the bands (at 95%) shown as dashed lines (green color), solid lines (red color) and the horizontal dotted (blue color).

5 Conclusion

From these simulation experiments and from the asymptotic theory, we draw the conclusion that the standard methodology, based on the QMLE, allows to fit SARMA representations of a wide class of nonlinear time series. But it is often restrictive to consider that the innovation process is directly observed. In future works, we intent to study how the existing estimation (see [16, 18]) and diagnostic checking (see [30]) procedures should be adapted in the situation where the GARCH process used in these simulation experiments is not directly observed, but constitutes the innovation of an observed SARMA-(seasonal)GARCH process which will be able to extend considerably the range of applications.

Appendix A Proofs

The proofs of Theorems 3 and 4 follow the same lines as in [7] and are similar. To have its own autonomy, the proofs will be rewrite and adapt.

Proof of Theorem 3

We recall that the Skorokhod space $\mathbb{D}^{k}[0{,}1]$ is the set of $\mathbb{R}^{k}-$valued functions defined on $[0{,}1]$ which are right continuous and have left limits. It is endowed with the Skorokhod topology and the weak convergence on $\mathbb{D}^{k}[0{,}1]$ is mentioned by $\xrightarrow[]{\mathbb{D}^{k}}$. We finally denote by $\lfloor x\rfloor$ the integer part of the real $x$.

To prove the result we need to recall some following results of [5, 13, 14, 15, 23].

We denote by $a_{i}^{*}$, $b_{i}^{*}$, $\mathbf{a}_{i}^{*}$ and $\mathbf{b}_{i}^{*}$ the coefficients defined by

[TABLE]

Following [13] and [23] (see also [5]), the noise derivatives involving in the expression of $J(\theta_{0})$ and $I(\theta_{0})$ can be represented as

[TABLE]

with $a_{i}^{*}=b_{i}^{*}=\mathbf{a}_{i}^{*}=\mathbf{b}_{i}^{*}=0$ when $i<0$.

For any $\theta\in\Theta\subset\mathbb{R}^{k_{0}}$ and any $(l,m)\in\{1,\dots,k_{0}\}^{2}$, under the above Assumptions, there exists absolutely summable and deterministic sequences $(c_{i}(\theta))_{i\geq 0}$, $(\lambda_{i,l}(\theta))_{i\geq 1}$ and $(\lambda_{i,l,m}(\theta))_{i\geq 1}$ such that, almost surely,

[TABLE]

with $c_{0}(\theta)=1$. A useful property of the above three sequences that they are asymptotically exponentially small. Indeed there exists $\rho\in]0{,}1[$ and a positive constant $K$ such that, for all $i\geq 1$, we have

[TABLE]

See Lemmas A.1. and A.2. of [15] for a more detailed treatment.

Now, in view of (12) it is clear that the asymptotic behaviour of $\hat{\gamma}_{m}$ is related to the limit distribution of $w_{t}=\left(-\Upsilon_{t}^{\prime}J^{-1^{\prime}}{,}\epsilon_{t}\epsilon_{t-1}{,}\dots,\epsilon_{t}\epsilon_{t-m}\right)^{\prime}$. First, we prove that $\frac{1}{\sqrt{n}}\sum_{j=1}^{\lfloor nr\rfloor}\Lambda w_{j}$ converges on the Skorokhod space to a Brownian motion. More precesily, we have to show that

[TABLE]

where $(B_{m}(r))_{r\geq 0}$ is a $m$-dimensional standard Brownian motion.

Using (19), the process $w_{t}$ can be rewritten as

[TABLE]

and thus the non-correlation between $\epsilon_{t}$ implies that $w_{t}$ has zero expectation with values in $\mathbb{R}^{k_{0}+m}$. In order to apply the functional central limit theorem for strongly mixing process, we need to identify the asymptotic covariance matrix in the classical central limit theorem for the sequence $(w_{t})_{t\geq 1}$. It is proved in Subsection 3.1 that

[TABLE]

where $f_{w}(0)$ is the spectral density of the stationary process $(w_{t})_{t\in\mathbb{Z}}$ evaluated at frequency 0. The main issue is to prove the existence of the matrix $\Xi$ which is a consequence (A3) and Davydov’s inequality [10]. For that sake, one has to introduce for any integer $k$, the random variables

[TABLE]

Since $w^{k}$ depends on a finite number of values of the noise-process $\epsilon$, it also satisfies a mixing property (see Theorem 14.1 in [9], p. 210). Based on the Davydov inequality (see [10]), the arguments developed in the Lemma A.1 in [13] (see also [14]) imply that

[TABLE]

where

[TABLE]

and thus (24) holds. Moreover we have that $\lim_{k\to\infty}\Xi_{k}=\Xi$.

Since the matrix $\Xi$ is positive definite, it can be factorized as $\Xi=\Delta\Delta^{\prime}$ where the $(k_{0}+m)\times(k_{0}+m)$ lower triangular matrix $\Delta$ has nonnegative diagonal entries. Therefore, we have

[TABLE]

and the new variance matrix can also been factorized as $\Lambda\Xi\Lambda^{\prime}=(\Lambda\Delta)(\Lambda\Delta)^{\prime}:=\Psi\Psi^{\prime}$, where $\Psi\in\mathbb{R}^{m\times(k_{0}+m)}$. Thus, ${n}^{-1/2}\sum_{t=1}^{n}\left(\Psi\Psi^{\prime}\right)^{-1/2}\Lambda w_{t}\xrightarrow[n\to\infty]{\mathrm{d}}\mathcal{N}\left(0,I_{m}\right)$ where $I_{m}$ is the identity matrix of order $m$. The above arguments also apply to matrix $\Xi_{k}$ with some matrix $\Psi_{k}$ which is defined analogously as $\Psi$. Consequently,

[TABLE]

and we also have ${n}^{-1/2}\sum_{t=1}^{n}\left(\Psi_{k}\Psi_{k}^{\prime}\right)^{-1/2}\Lambda w^{k}_{t}\xrightarrow[n\to\infty]{\mathrm{d}}\mathcal{N}\left(0,I_{m}\right)$

Now we are able to apply the functional central limit theorem for strongly mixing process of [19]. We have for any $r\in(0{,}1)$,

[TABLE]

For all $j\in\{1,\dots,[nr]\}$, we write

[TABLE]

and we obtain that

[TABLE]

In order to conclude (23), it remains to observe that, uniformly with respect to $n$,

[TABLE]

where

[TABLE]

By Lemma 4 in [14], we have

[TABLE]

and since $\lfloor nr\rfloor\leq n$,

[TABLE]

Thus (26) is true and the proof of the first step (23) is achieved.

The previous step ensures us that Assumption 1 in [22] is satisfied for the sequence $(\Lambda w_{t})_{t\geq 1}$. We follow the arguments developed in Sections 2 and 3 in [22], the second step is to show that

[TABLE]

by applying the continuous mapping theorem on the Skorokhod space and where the random variable $V_{m}$ is defined in (13). The main issue is to obtain that

[TABLE]

by continuous mapping theorem and using (23), the fact that $[nr]/n\to r$ as $n\to\infty$. In view of (28), it follows that

[TABLE]

which prove (27). Since $\sqrt{n}\hat{\gamma}_{m}=n^{-1/2}\sum_{t=1}^{n}\Lambda w_{t}+\mathrm{o}_{\mathbb{P}}(1)$, using (23) and (27) we obtain

[TABLE]

The proof of Theorem 3 is then complete.     $\Box$

Proof of Theorem 4:

We write $\hat{C}_{m}=C_{m}+\nabla_{n}$ where \nabla_{n}={n}^{-2}\sum_{t=1}^{n}\big{(}S_{t}S_{t}^{\prime}-\hat{S}_{t}\hat{S}_{t}^{\prime}\big{)}. There are three kinds of entries in the matrix $\nabla_{n}$. The first one is a sum composed of

[TABLE]

for $(k,k^{\prime})\in\{1,\dots,m\}^{2}$. Using (22) and the consistency of $\hat{\theta}_{n}$, we have $\upsilon_{t}^{k,k^{\prime}}=\mathrm{o}(1)$ almost surely. The two last kinds of entries of $\nabla_{n}$ come from the following quantities for $i,j\in\{1,\dots,k_{0}\}$ and $k\in\{1,\dots,m\}$

[TABLE]

and they also satisfy $\tilde{\upsilon}_{t}^{k,i}+\bar{\upsilon}_{t}^{i,j}=\mathrm{o}(1)$ almost surely by using (19) and (22). Consequently, $\nabla_{n}=\mathrm{o}(1)$ almost surely as $n$ goes to infinity. Thus one may find a matrix $\nabla^{\ast}_{n}$, that tends to the null matrix almost surely, such that

[TABLE]

Thanks to the arguments developed in the proof of Theorem 3, $n\hat{\gamma}_{m}^{\prime}C_{m}^{-1}\hat{\gamma}_{m}$ converges in distribution. So $n\hat{\gamma}_{m}^{\prime}\nabla_{n}^{\ast}\hat{\gamma}_{m}$ tends to zero in distribution, hence in probability. Then $n\hat{\gamma}_{m}^{\prime}\hat{C}_{m}^{-1}\hat{\gamma}_{m}$ and $n\hat{\gamma}_{m}^{\prime}C_{m}^{-1}\hat{\gamma}_{m}$ have the same limit in distribution and the result is proved.      $\Box$

Acknowledgements

We sincerely thank the anonymous reviewers and Editor for helpful remarks. The authors wish to acknowledge the support from the "Séries temporelles et valeurs extrêmes : théorie et applications en modélisation et estimation des risques" Projet Région (Bourgogne Franche-Comté, France) grant No OPE-2017-0068.