State space time series modelling of the Dutch Labour Force Survey: Model selection and mean squared errors estimation
Section 2. The Dutch Labour Force Survey

• Previous
• Next

2.1 The DLFS design

The DLFS has been based on a rotating panel design since October 1999. Every month, a sample of addresses is drawn according to a stratified two-stage sample design. Strata are formed by geographical regions; municipalities are the primary sampling units and addresses are the secondary sampling units. All households residing on one address are included in the sample. In this paper, the DLFS data observed from January 2001 until June 2010 are considered. During this period, data in the first wave were collected by means of computer assisted personal interviewing (CAPI) by interviewers that visit sampled households at home. After a maximum of six attempts, an interviewer leaves a letter with the request to contact the interviewer by telephone to make an appointment for an interview. When a household member cannot be contacted, proxy interviewing is allowed by members of the same household. Respondents are re-interviewed four times at quarterly intervals. In these four subsequent waves, data are collected by means of computer assisted telephone interviewing (CATI). During these re-interviews, a condensed questionnaire is applied to establish any changes in the labour market position of the respondents. Proxy interviewing is also allowed during these re-interviews. Mobile phone numbers and secret land line numbers are collected in the first wave to avoid panel attrition. With the commencement of the rotating panel design for the DLFS, the gross sample size was about 6,200 addresses per month on average, with about 65% completely responding households. The response rates in the follow-up waves are about 90% compared to the preceding wave.

The general regression (GREG) estimator (Särndal et al. 1992) is applied to estimate the total unemployed labour force. This estimator accounts for the complexity of the sample design and uses auxiliary information available from registers to correct, at least partially, for selective non-response. Let $Y_t^j$ denote the GREG estimate of the total number of unemployed in month $t$ based on the $j^{\text{th}}$ wave of respondents. Five such estimates are obtained per month, each of them being respectively based on the sample that entered the survey in month $t-l,$ $l\mathrm{<mo>=</mo>}\left\{\mathrm{0,}\mathrm{3,}\mathrm{6,}\mathrm{9,}12\right\}.$ The GREG estimator for this population total is defined as:

$$Y_t^j\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\in s}{w}_{k\mathrm{,}t}}\left({\displaystyle \sum _{i=1}^{n_{k,t}}{y}_{i\mathrm{,}k\mathrm{,}t}}\right)(2.1)$$

with $y_{i\mathrm{,}k\mathrm{,}t}$ representing the sample observations that are equal to 1 if the $i^{\text{th}}$ person in the $k^{\text{th}}$ household is unemployed, and zero otherwise; $n_{k\mathrm{,}t}$ is the number of persons aged 15 or above in the $k^{\text{th}}$ household; $w_{k\mathrm{,}t}$ are the regression weights for household $k$ at time $t.$ The method of Lemaître and Dufour (1987) is used to obtain equal weights for all persons within the same household:

$$w_{k\mathrm{,}t}\mathrm{<mo>=</mo>}\frac{1}{{\pi }_{k\mathrm{,}t}}\left[1+\left(X_t-{\displaystyle \sum _{k\in s}}\frac{x_{k,t}}{{\pi }_{k\mathrm{,}t}}\right){\left({\displaystyle \sum _{k\in s}}\frac{x_{k,t}{x}_{k,t}^{{}^{\prime }}}{{\pi }_{k\mathrm{,}t}{g}_{k\mathrm{,}t}}\right)}^{-1}\frac{x_{k,t}}{g_{k\mathrm{,}t}}\right]\mathrm{,}(2.2)$$

where ${\pi }_{k\mathrm{,}t}$ is the inclusion probability of household $k$ at time $t,$ $g_{k\mathrm{,}t}$ is the size of household $k$ at time $t;$ $x_{k,t}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i=1}^{n_{k,t}}{x}_{i,k,t}},$ with $x_{i,k,t}$ being a $J-$ dimensional vector with the weighting model auxiliary information on the $i^{\text{th}}$ person in the $k^{\text{th}}$ household at time $t.$ Vector $X_t$ contains population totals of auxiliary variables. The weighting model is defined by the following variables (with the number of categories in brackets): Age(5)Gender + Geographic Region(44) + Gender(2) $\times$ Age(21) + Age(5) $\times$ Marital Status(2) + Ethnicity(8), where $\times$ stands for interaction of variables, and Age(5)Gender is a variable classified into eight classes where Age has five categories, with the second, third and fourth age category being itemized into two genders.

The variance of the GREG estimator $Y_t^j$ is approximated by:

$$\widehat{\text{Var}}\left(Y_t^j\right)\mathrm{<mo>=</mo>}{\displaystyle \sum _{h\mathrm{<mo>=</mo>\; 1}}^H}\frac{n_{h\mathrm{,}t}}{n_{h\mathrm{,}t}-1}\left({\displaystyle \sum _{k=1}^{n_{h,t}}{\left(w_{k\mathrm{,}t}{\widehat{e}}_{k\mathrm{,}t}\right)}^2}-\frac{1}{n_{h\mathrm{,}t}}{\left({\displaystyle \sum _{k=1}^{n_{h,t}}{w}_{k\mathrm{,}t}{\widehat{e}}_{k\mathrm{,}t}}\right)}^2\right)\mathrm{,}j\mathrm{<mo>=</mo>\; <mo>\{</mo>\; 1,}\mathrm{2,}\mathrm{3,}\mathrm{4,}\mathrm{5\; <mo>\}</mo>\; ,}(2.3)$$

where the GREG residuals are ${\widehat{e}}_{k\mathrm{,}t}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i=1}^{n_{k,t}}\left(y_{i\mathrm{,}k\mathrm{,}t}-x_{i,k,t}^{{}^{\prime }}{\widehat{\beta }}_t\right)};$ $n_{h\mathrm{,}t}$ is the number of households in stratum $h$ (with $H$ being the total number of strata); vector ${\widehat{\beta }}_t$ is a Horvitz-Thompson type estimator for the regression coefficient that is obtained from regressing the target variable on the auxiliary variables from the sample.

2.2 The STS model for the DLFS

There are two reasons why Statistics Netherlands took a decision to switch to a time series model-based production approach in June 2010. One reason for that was inadequately small sample sizes for production of monthly estimates. With a net sample size of about 4,000 households in the first wave on average, the GREG estimates of the unemployed labour force had a coefficient of variation of about 4% at the national level, which was considered to be too volatile for official statistical publications. In addition to that, monthly unemployment figures must be published for six domains based on a classification of gender and age. The design-based estimates of these domains feature much higher coefficients of variation. Another problem with the DLFS is the so-called RGB, which refers to systematic differences between the estimates of different waves (see, e.g., Bailar 1975 or Pfeffermann 1991). Common reasons behind the RGB are panel attrition, panel-effects, and differences in questionnaires and modes used in the subsequent waves. In the case of the DLFS, the first wave estimates are assumed to be most reliable, with the subsequent waves systematically underestimating the unemployed labour force numbers. See van den Brakel and Krieg (2009) for a more detailed discussion.

Both problems are solved with an STS model, which uses five series of GREG estimates for the five different waves as input. With an STS model, an observed series is decomposed into several unobserved components, e.g., trend and seasonal. The Kalman filter, optionally in combination with a smoothing algorithm, can be applied to extract these components from the observed time series. By doing so, estimates of the components that define the signal for unemployment are separated from unexplained variance of the population parameter and from the sampling variance. This generally results in less volatile point estimates, with substantially smaller standard errors compared to those of the GREG estimates. By modelling the systematic differences between the five input series, the model also accounts for the RGB of the rotating panel.

In each month $t,$ a five-dimensional vector $Y_t\mathrm{<mo>=</mo>}{\left(Y_t^1{Y}_t^2{Y}_t^3{Y}_t^4{Y}_t^5\right)}^{\prime }$ is observed, containing GREG estimates of the total number of the unemployed labour force based on the five waves. Based on Pfeffermann (1991), van den Brakel and Krieg (2009) developed the following model for the GREG estimates $Y_t:$

$$Y_t=1_5{\xi }_t+{\lambda }_t+e_t\mathrm{,}(2.4)$$

here, $1_5$ is a five-dimensional column vector of ones, ${\xi }_t$ is the unknown (scalar) true population parameter, ${\lambda }_t$ is a vector containing state variables for the RGB, and $e_t$ is a vector of the survey errors that are correlated with their counterparts from previous waves (the structure will be presented later). For the true population parameter, it is assumed that: ${\xi }_t\mathrm{<mo>=</mo>}{L}_t+{\gamma }_t+{\epsilon }_t,$ which is the sum of a stochastic trend $L_t,$ a stochastic seasonal component ${\gamma }_t,$ and an irregular component ${\epsilon }_t\stackrel{\text{iid}}{\mathrm{~}}N\left(\mathrm{0,}{\sigma }_{\epsilon }^2\right).$

For the stochastic trend $L_t,$ the so-called smooth-trend model is assumed:

$$\begin{array}{ll}{L}_t\hfill & \mathrm{<mo>=</mo>}{L}_{t-1}+R_{t-1}\mathrm{,}\hfill \\ R_t\hfill & \mathrm{<mo>=</mo>}{R}_{t-1}+{\eta }_{R\mathrm{,}t}\mathrm{,}\hfill \end{array}$$

where $L_t$ and $R_t$ represent the level and slope of the true population parameter, respectively, with the slope disturbance term being distributed as: ${\eta }_{R\mathrm{,}t}\stackrel{\text{iid}}{\mathrm{~}}N\left(\mathrm{0,}{\sigma }_R^2\right).$

For the seasonal component ${\gamma }_t,$ the trigonometric model is assumed:

$${\gamma }_t\mathrm{<mo>=</mo>}{\displaystyle \sum _{l\mathrm{<mo>=</mo>\; 1}}^6{\gamma }_{t\mathrm{,}l}}\mathrm{,}$$

where each of these six harmonics follows the process:

$$\begin{array}{ll}{\gamma }_{t\text{}\mathrm{,}l}\hfill & \mathrm{<mo>=</mo>}\text{cos}\left(h_l\right){\gamma }_{t-\mathrm{1,}l}+\text{sin}\left(h_l\right){\gamma }_{t-\mathrm{1,}l}^{*}+{\omega }_{t\text{}\mathrm{,}l}\mathrm{,}\hfill \\ {\gamma }_{t\text{}\mathrm{,}l}^{*}\hfill & \mathrm{<mo>=</mo>}-\text{sin}\left(h_l\right){\gamma }_{t-\mathrm{1,}l}+\text{cos}\left(h_l\right){\gamma }_{t-\mathrm{1,}l}^{*}+{\omega }_{t\text{}\mathrm{,}l}^{*}\mathrm{,}\hfill \end{array}$$

with $h_l\mathrm{<mo>=</mo>}\frac{\pi l}{6}$ being the $l^{\text{th}}$ seasonal frequency, $l\mathrm{<mo>=</mo>}\left\{1,\dots 6\right\}.$ The zero-expectation stochastic terms ${\omega }_{t\mathrm{,}l}$ and ${\omega }_{t\mathrm{,}l}^{*}$ are assumed to be normally and independently distributed and to possess the same variance within and across all the harmonics, such that:

$$\begin{array}{ll}\text{Cov}\left({\omega }_{t\mathrm{,}l}\mathrm{,}{\omega }_{t^{\prime },l^{\prime }}\right)\hfill & \mathrm{<mo>=</mo>}\text{Cov}\left({\omega }_{t\mathrm{,}l}^{*}\mathrm{,}{\omega }_{t^{\prime },l^{\prime }}^{*}\right)\mathrm{<mo>=</mo>}(\begin{array}{ll}{\sigma }_{\omega }^2\hfill & \text{if}l=l^{\prime }\text{and}t=t^{\prime },\hfill \\ 0\hfill & \text{if}l\ne l^{\prime }\text{or}t\ne t^{\prime },\hfill \end{array}\hfill \\ \text{Cov}\left({\omega }_{t\mathrm{,}l}\mathrm{,}{\omega }_{t\mathrm{,}l}^{*}\right)\hfill & \mathrm{<mo>=</mo>\; 0}\text{for}\text{all}l\text{and}t\mathrm{.}\hfill \end{array}$$

The second component in (2.4) is the RGB. It is assumed that the first wave is unbiased, as motivated in van den Brakel and Krieg (2009). The RGBs for the follow-up waves are time-dependent and are modelled as random walk processes. The rationale behind this is that field-work procedures are subject to frequent changes. Apart from that, response rates change gradually over time. This makes the RGB time-dependent, as illustrated by van den Brakel and Krieg (2015), Figure 4.3. The RGB vector for the five waves can be written in the following form: ${\lambda }_t\mathrm{<mo>=</mo>}{\left(0{\lambda }_t^2{\lambda }_t^3{\lambda }_t^4{\lambda }_t^5\right)}^{\prime },$ with:

$${\lambda }_t^j\mathrm{<mo>=</mo>}{\lambda }_{t-1}^j+{\eta }_{\lambda \mathrm{,}t}^j\mathrm{,}j\mathrm{<mo>=</mo>}\left\{\mathrm{2,}\mathrm{3,}\mathrm{4,}5\right\}\mathrm{.}$$

It is assumed that the RGB disturbances are not correlated across different waves and are normally distributed: ${\eta }_{\lambda \mathrm{,}t}^j\stackrel{\text{iid}}{\mathrm{~}}\left(\mathrm{0,}{\sigma }_{\lambda }^2\right),$ with equal variances in all the four waves.

The last component in (2.4) contains the survey errors for the five GREG estimates, i.e., $e_t\mathrm{<mo>=</mo>}{\left(e_t^1{e}_t^2{e}_t^3{e}_t^4{e}_t^5\right)}^{\prime }.$ To account for sampling error heterogeneity caused by changes in the sample sizes over time, the sampling errors are modelled proportionally to the design-based standard errors according to the following measurement error model proposed by Binder and Dick (1990): $e_t^j\mathrm{<mo>=</mo>}{\tilde{e}}_t^j{z}_t^j,$ where $z_t^j\mathrm{<mo>=</mo>}\sqrt{\widehat{\text{Var}}\left(Y_t^j\right)}$ and ${\tilde{e}}_t^j$ are standardised sampling errors that follow a stationary process defined later in the text. Here, $\widehat{\text{Var}}\left(Y_t^j\right)$ are the design-based variance estimates obtained from the micro data using (2.3). They are treated as a priori known sampling variances in the STS model.

Since the sample in the first wave has no overlap with samples observed in the past, ${\tilde{e}}_t^t$ can be modelled as a white noise with $E\left({\tilde{e}}_t^1\right)\mathrm{<mo>=</mo>\; 0}$ and $\text{Var}\left({\tilde{e}}_t^1\right)\mathrm{<mo>=</mo>}{\sigma }_{v_1}^2.$ The variance of the survey errors $e_t^t$ will be equal to the variance of the GREG estimates if the maximum likelihood estimate of ${\sigma }_{v_1}^2$ is approximately equal to unity.

The survey errors in the follow-up waves are correlated with the survey errors from the preceding waves. This autocorrelation coefficient is estimated from the survey data using the approach proposed by Pfeffermann, Feder and Signorelli (1998). The autocorrelation structure is modelled with an AR(1) model where the autocorrelation coefficient is obtained with the Yule-Walker equations (van den Brakel and Krieg 2009):

$${\tilde{e}}_t^j\mathrm{<mo>=</mo>}\rho {\tilde{e}}_{t-3}^{j-1}+{\nu }_t^j\mathrm{,}{\nu }_t^j\stackrel{\text{iid}}{\mathrm{~}}N\left(\mathrm{0,}{\sigma }_{v_j}^2\right)\mathrm{,}j\mathrm{<mo>=</mo>}\left\{\mathrm{2,<mtext>\hspace{0.17em}</mtext>3,<mtext>\hspace{0.17em}</mtext>4,<mtext>\hspace{0.17em}</mtext>5}\right\}\mathrm{.}$$

It is assumed that the first-order autocorrelation coefficient is common for all the four waves. Its estimate is used as a priori information in the model. Since ${\tilde{e}}_t^j$ is an AR(1) process, $\text{Var}\left({\tilde{e}}_t^j\right)\mathrm{<mo>=</mo>}{\sigma }_{v_j}^2/\left(1-{\rho }^2\right).$ The variance of the sampling error $e_t^j$ is approximately equal to $\widehat{\text{Var}}\left(Y_t^j\right)$ if the maximum likelihood estimate of ${\sigma }_{v_j}^2$ is approximately equal to $\left(1-{\rho }^2\right).$ Five different hyperparameters ${\sigma }_{v_j}^2\text{}\mathrm{,}j\mathrm{<mo>=</mo>}\left\{\mathrm{1,<mtext>\hspace{0.17em}</mtext>2,<mtext>\hspace{0.17em}</mtext>3,<mtext>\hspace{0.17em}</mtext>4,<mtext>\hspace{0.17em}</mtext>5}\right\},$ are assumed for the survey error components of the five waves.

The disturbance variances, together with the autocorrelation parameter $\rho ,$ are collected in a hyperparameter vector called $\theta \mathrm{<mo>=</mo>}{\left({\sigma }_R^2{\sigma }_{\omega }^2{\sigma }_{\epsilon }^2{\sigma }_{\lambda }^2{\sigma }_{v_1}^2{\sigma }_{v_2}^2{\sigma }_{v_3}^2{\sigma }_{v_4}^2{\sigma }_{v_5}^2\rho \right)}^{\prime },$ and the vector containing only the disturbance variances is called ${\theta }_{\sigma }\mathrm{<mo>=</mo>}{\left({\sigma }_R^2{\sigma }_{\omega }^2{\sigma }_{\epsilon }^2{\sigma }_{\lambda }^2{\sigma }_{v_1}^2{\sigma }_{v_2}^2{\sigma }_{v_3}^2{\sigma }_{v_4}^2{\sigma }_{v_5}^2\right)}^{\prime }.$ To avoid negative estimates, the disturbance variance hyperparameters in ${\theta }_{\sigma }$ are estimated on a log-scale. The quasi-maximum likelihood method is used (see e.g., Harvey 1989), where $\widehat{\rho }-$ estimates are treated as known. Numerical analysis of this paper is conducted with OxMetrics 5 (Doornik 2007) in combination with SsfPack 3.0 package (Koopman, Shephard and Doornik 2008).

• Previous
• Article main page
• Next