Model-Based Seasonal Adjustment Diagnostics

Tucker McElroy

U.S. Census Bureau

Abstract

Several model-based seasonal adjustment diagnostics are currently being incorporated into X-13-ARIMASEATS, a new hybrid seasonal adjustment program under development at the U.S. Census Bureau. These

diagnostics have been implemented in the Ox programming language, and this paper discusses various

aspects of the implementation. The diagnostics are normalized quadratic forms in the data, and roughly

correspond to intra- and inter-component variation.

Disclaimer This report is released to inform interested parties of research and to encourage discussion.

The views expressed on statistical issues are those of the author and not necessarily those of the U.S. Census

Bureau.

1 Introduction

A natural way to assess the quality of signal extraction is to check the spectrum of the estimated noise

for residual signal. However, it is difficult to quantify any given form of assessment without making further

assumptions about the process. One approach is to stipulate a model for the data, and gauge the degree to

which the spectrum for the estimated noise deviates from what the model predicts. In the model-based signal

extraction scenario such an approach – formulated in the time domain – is presented in Maravall (2003),

which first proposed a single statistic that would quantify the goodness of signal extraction. Maravall’s

diagnostic was earlier incorporated into SEATS, the widely-used model-based seasonal adjustment program

of G ́omez and Maravall (1994), but wasn’t documented until Maravall (2003). The basic concept of the diagnostic was modified and studied in Findley, McElroy, and Wills (2004), adapting the idea for finite sample

signal extraction. The concept of the intra-component diagnostic is to measure the variation of an estimated

component, assessed through a variance estimate of the appropriately “differenced” signal extraction estimate, and compare this quantity to what we would expect if our model were true. The inter-component

diagnostic assesses variability between estimated components in the same manner. Extreme values of the

statistic indicate model inadequacy – and hence filter inadequacy – with respect to the component model

for the desired signal. Findley et al. (2004) discuss the basic theory of the new signal extraction diagnostics,

and describe the results of simulation and empirical studies, while McElroy (2005b) expands and generalizes

these ideas, providing asymptotic results. The paper at hand relates the implementation of these methods in

the Ox (Doornik, 1998) program SigDiagObj.ox, utilizing the SsfPack function suite (Koopman, Shepherd,

and Doornik, 1999).

1
Findley, et al. (2004) principally dealt with simple seasonal models, and focused on the irregular component. The current implementation considers a greater variety of components, and a broader class of models.

In particular, we consider data decompositions of the form

Yt = St + Tt + It, (1)

where the data is Yt, the seasonal is St, the trend-cycle is Tt, and the irregular is It. The program SigDiagObj.ox conducts four basic operations: estimation, decomposition, signal extraction, and diagnostic

calculation. This paper discusses each of these steps in some detail, illustrated through the series of U.S.

Retail Sales of Shoe Stores data from the monthly Retail Trade Survey of the Census Bureau, from 1984 to

1998, which will be referred to as the “shoe” series.

2 Estimation

In order to do model-based signal extraction, it is necessary to specify models for each of the components.

We will follow the canonical decomposition approach of Hillmer and Tiao (1982), and so we start with a

Seasonal ARIMA (SARIMA) model for monthly data Yt, after suitable Box-Cox transformations:

φ(B)Φ(B

12)(1 − B)

d

(1 − B

12)

D

Yt = θ(B)Θ(B

12)2Y

t

t = 1, 2, · · · ,n (2)

where φ(z), Φ(z), θ(z), and Θ(z) are polynomials with roots outside the unit circle of the complex plane,

ensuring an invertible representation for the differenced data. Here 2

Y

t

is a white noise innovation sequence.

We allow for d nonseasonal differences and D seasonal differences. Currently the program works with d, D ≤ 1

but not both zero; this is not a real restriction in practice, since these differencing orders allow for a fairly

wide range of nonstationary processes. The assumed component models are given as follows:

Φ(B

12)U(B)

D

St = θ

S

(B)2t

S

(3)

φ(B)(1 − B)

d+D

Tt = θ

T

(B)2

T

t

It = 2t

I

,

where the various 2t sequences are independent white noise. As above, the MA polynomials θ

S(z) and θT

(z)

are chosen so as to guarantee an invertible representation. The nonseasonal AR polynomial is associated

with the trend-cycle component Tt, and the seasonal AR goes with the seasonal component. The polynomial

U(z) = 1+z +z

2 +· · ·+z11 achieves seasonal differencing, and satisfies 1−z12 = U(z)(1−z). The irregular

component is assumed to be white noise. The canonical decomposition method begins by estimating all the

parameters (including the innovation variance σ

2

) of the data model (??).

2Y The most popular estimates

are computed by the method of maximum likelihood, whose implementation in SsfPack implicitly assumes a

joint Gaussian distribution for the data. SigDiagObj.ox uses SsfPack’s estimation method, which maps the

SARIMA model (??) into State Space Form, and maximizes the likelihood using state space methods. This

is numerically efficient and stable, and is to be preferred over more direct approaches. One limitation is that

SigDiagObj.ox will not select model orders for the user – one must know ahead of time what SARIMA model

orders are desired. Essentially, the user specifies the differencing orders (d, D) and a collection of initial

parameter values, which implicitly determine the SARIMA model orders (p, P, q, Q) for the autoregressive

and moving average polynomials.

2
One is also able to fix all of the parameters with user-selected values. This facility is useful for testing a

priori defined filters (i.e., model-based filters whose parameters are not determined by the data, but rather by

the user) on the data. One snag of the Hillmer-Tiao method is that some SARIMA models are inadmissable

in the sense that they do not possess a decomposition. By overwriting certain parameter values, the user

may be able to impose an admissible model on the data.

Example Suppose we specify a SARIMA(1, 1, 1)(0, 1, 1)12 model for the logged shoe series. We obtain

the following data model via estimation:

(1 − .15B)(1 − B)(1 − B

12)Yt = (1 − .67B)(1 − .35B12)2Y

, σ2

t = .00095 2Y

This is accomplished through calling InitParams (< .6 >, <>, < −.6 >, < −.5 >, 0,< 1, 1 >) – which specifies

the model and initializes the parameters – and Estimate(), which uses maximum likelihood estimation to

produce parameter estimates. More detailed information on InitParams is provided in its documentation

within SigDiagObj.ox.

3 Decomposition

We first stipulate component models according to (??), a process that involves assigning all the “lefthand” operators of (??) to the left-hand sides of the component models. These left-hand operators include

the AR operators and the differencing operators; for the latter, it is important that the factors are distributed

uniquely, so that the differencing polynomials for each component are relatively prime. This is a prerequisite

of the model-based signal extraction theory – see McElroy (2005a). The MA operators for the component

models are then determined according to the canonical decomposition method of Hillmer and Tiao (1982),

implemented through the Ox routine gendecomp.ox of Aston and Koopman (2003).

The gendecomp.ox function operates in the same manner as the code in SEATS for producing canonical

decompositions, although our implementation is somewhat more simplistic. By associating all nonseasonal

AR operators with the trend, we essentially define a trend-cycle component, which is more simply referred

to as the trend. In contrast, SEATS will often generate a fourth cycle component that is an ARMA process,

with certain of the nonseasonal AR factors associated to it. We have chosen to lump the trend and cycle

together, which seems to be a reasonable approach given our principal objective of assessing the quality of

seasonal adjustments. If one is further interested in cyclical behavior, a four-component decomposition could

be attempted, though this may encounter problems of admissability.

Example In the shoe series, the AR(1) operator is associated with the trend:

U(B)St = (1 + 1.26B + 1.06B

2 + .84B3 + .58B4 + .35B5 + .15B6

− .03B

7 − .14B8 − .26B9 − .25B10 − .42B11)2S

t

, σ2

2S

(1 − .15B)(1 − B)

2

Tt = (1 − .46B − .96B

2 + .49B3

)2

T

t

, σ2

2

T = .000053

= .00012

It = 2t

I

, σ2

2

I = .00023

The above results were produced by calling the Decompose() function; the PrintComponents() function may

be used to output the models to the screen.

3
4 Signal Extraction

We adopt the model-based, finite sample approach to signal extraction. The various estimates can be

produced through a state space smoother, implemented for example in SsfPack, but we produce the filter

matrix directly utilizing formulas in McElroy (2005a). This is necessary for our subsequent calculations of

the diagnostics themselves – see Findley et al. (2004). For each of the three components – trend, seasonal,

and irregular – we produce the appropriate signal extraction matrix, each of whose rows corresponds to the

time-dependent filter for the corresponding time point. In other words, the ith row of the filter matrix F

consists of the filter coefficients that, when applied to the data Yt, produce an estimate of the signal at time

i. For example, if estimating the seasonal we have

FY = Sˆ = E[S|Y ],

0 0

where Y = (Y1,Y2, · · · ,Yn) and S = (S1, S2, · · · ,Sn) is the seasonal. In this stage of the program, we

produce the filter matrix F and the signal estimate Sˆ.

There are six quantities considered as possible signals of interest. These represent appropriately differenced

components or combinations of such as follows:

U(B)

D

St, (1 − B)

d+D

Tt, It

U(B)

D

(St + It), (1 − B)

d+D

(Tt + It), (1 − B)

d

(1 − B

12)

D

(St + Tt)

The basic intra-component diagnostics are sample autocovariances – computed at various lags – of the differenced signal estimates, which means that we first estimate the differenced signals U(B)

D

St, (1 − B)

d+D

Tt,

etc., and then compute the lagged average of squares. If u denotes a vector of estimates for a differenced

signal, then u

0u/n gives the lag zero sample autocovariance, assuming mean zero. If we let L be a lag matrix

(of dimension equal to the length of u) with Lij equal to 1 if i = j + 1 and 0 otherwise, then u

0Lhu/n is a

quadratic form that yields the lag h estimate of the autocovariance (h is a non-negative integer).

The inter-component diagnostics are sample cross-covariances, computed for six pairings of components:

S, I; S, T; S, TI; T,I; T, SI; and TS, I. We have a similar calculation, but now it involves two vectors u and

v, containing estimates for two diverse differenced signals. An additional complication is that u and v could

have different lengths (e.g., the irregular has length n but differenced trend has length n − (d + D)). Based

on theoretical considerations, it is appropriate to trim the longer vector of its first values such that the two

vectors have the same length; call this v ̃ if v was the longer vector. Then u

0v ̃/n gives the cross-covariance

estimate.

Example In order to compute the signal extraction filters, we first call buildDiffMatrices() and then

buildCovMatrices(), which construct the differencing and covariance matrices that will be needed. Then

ExtractSignals() computes the filters and applies them to the data. Below in Figure 1 is a picture of the

trend-cycle estimate in the shoe series together with the data.

4
Figure 1: Logarithm of U.S. Retail Sales of Shoe Stores, with Trend-Cycle estimate

Year

1984 1986 1988 1990 1992 1994 1996 1998

6.8 7.0 7.2 7.4 7.6

Data

Trend−Cycle

5 Diagnostic Calculation

The diagnostics mentioned above provide a measure of variation in an estimated signal. Under the Null

Hypothesis that our model for Yt is correct, we can compute the mean and variance of this diagnostic and

thereby construct a standardized diagnostic. Further details are developed in McElroy (2005b); here we

discuss the interpretation and application of the diagnostics.

The standardized diagnostic is asymptotically normal under some mild conditions, and hence large positive

or negative values indicate rejection of the specified model, whether this was generated through estimation

or a priori methods. Since some of the diagnostics may be significant while others are not, the extreme

values indicate model inadequacy only in certain portions of the data’s spectrum. Roughly speaking, a

significant diagnostic indicates that the spectral density of the differenced data Wt = (1 − B)

d

(1 −

D B12) Yt

is poorly modelled at the frequencies corresponding to the component of interest. For example, an extreme

value of the trend diagnostic indicates poor modelling (and hence filtering) of the low frequencies, whereas

an extreme value of the seasonal diagnostic would show poor modelling of the six seasonal frequencies. This

interpretation is not completely rigorous – the diagnostics need to be further modified for that to happen.

It is also difficult to interpret positive versus negative values of the diagnostic; it is safer to stick to the

two-sided alternative. See McElroy (2005b) for a discussion.

The final output of the program includes the p-values for a two-sided test along with the standardized

diagnostics, for each of the six signals defined above. Significant p-values indicate that adjustment to the

model may be needed to improve signal extraction.

5
Example We call the function ComputeLagDiagnostics(h) with h = 0, 1, 12 and obtain the following

output:

Table 1. Auto-covariance Diagnostics with p-values

Signal Lag Zero p-value Lag One p-value Lag Twelve p-value

Irregular −0.423 .336 0.226 .411 0.648 .258

Trend 0.088 .465 −1.168 .121 0.631 .264

Seasonal −0.043 .483 −0.058 .477 0.214 .415

Trend-Irregular −0.467 .320 0.697 .243 0.458 .324

Seasonal-Irregular −0.082 .467 0.019 .492 0.101 .460

Trend-Seasonal −0.213 .416 0.450 .326 −0.171 .432

These diagnostics indicate adequacy of the filters, since none of the p-values are significant. For crosscovariance, we call the function ComputeCrossDiagnostics() and obtain:

Table 2. Cross-covariance Diagnostics with p-values

Signal Pair Diagnostic p-value

Seasonal, Irregular 0.637 .262

Seasonal, Trend −0.853 .197

Seasonal, Trend-Irregular −1.412 .079

Trend, Irregular 1.067 .143

Trend, Seasonal-Irregular −0.345 .365

Trend-Seasonal, Irregular 0.738 .230

Again, the high p-values indicate model adequacy with respect to filtration.

No-Ends In Findley, McElroy, and Wills (2004), a further adjustment to the diagnostic is suggested,

namely to trim the beginning and ends of each vector of estimates u, and adjust the statistical normalization

appropriately. In that paper, it was proposed that this “no-ends” version would eliminate irritating end

effects that interfered with the finite sample performance of the diagnostic. SigDiagObj.ox can also produce

the “no-ends” version of the above diagnostics (not shown here).

6 Conclusion

The program SigDiagObj.ox performs SARIMA model estimation, canonical decomposition, signal extraction, and diagnostic calculation, given a model specification. Current implementation assumes a decomposition into seasonal, trend-cycle, and irregular; we hope to generalize this to four components if practicable.

Also, the constraints on the differencing orders d and D make the current program easier to apply; future

work will focus on addressing these constraints in the interest of making the program more broadly applicable. The incorporation of these diagnostics into X–13–ARIMA–SEATS will allow their testing upon

hundreds of series at once, facilitating a large-scale empirical study. Also under development is a version

of SigDiagObj.ox that allows for a separate “transient” component, which can model a separate cycle or a

sampling error component, for example. This new version, called Hybrid.ox, was designed to be an extension

6
of the diagnostics to handle series with sampling error, such as those encountered at the Bureau of Labor

Statistics.

References

[1] Aston, J. and Koopman, S. (2003) An Implementation of Component Models for Seasonal Adjustment

Using the SsfPack Software Module of Ox. Proceedings of the 13th Federal Forecasters Conference, 47–53.

[2] Doornik, J. (1998) Object-oriented Matrix Programming using Ox 2.0. London: Timberlake Consultants

Press.

[3] Findley, D., McElroy, T., and Wills, K. (2004) Modification of SEATS’ Diagnostic for Detecting Overor Underestimation of Seasonal Adjustment Decomposition Components. 2004 Proceedings American

Statistical Association [CD-ROM]: Alexandria, VA.

[4] G ́omez, V. and Maravall, A. (1994) Program SEATS Signal Extraction in ARIMA Time Series: Instructions for the User. Working Paper ECO 94/28, European University Institute, Florence.

[5] Hilmer, S. and Tiao, G. (1982) An ARIMA-model-based approach to seasonal adjustment. Journal of

the American Statistical Association 77, 377, 63–70.

[6] Koopman, S., Shepherd, N., and Doornik, J. (1999) Statistical Algorithms for Models in State Space

Using SsfPack 2.2. Econometrics Journal 2, 113–166.

[7] Maravall, A. (2003) A Class of Diagnostics in the ARIMA-Model-Based Decomposition of a Time Series.

Memorandum, Bank of Spain.

[8] McElroy, T. (2005a) Matrix Formulas for Nonstationary Signal Extraction. SRD Research Report No.

RRS2005 − 04, Bureau of the Census.

[9] McElroy, T. (2005b) Statistical Properties of Signal Extraction Diagnostics. Working Paper.

7