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

• Previous
• Next

5.1 Alternative model specifications for the DLFS

STS models are usually selected and evaluated by means of formal diagnostic tests for normality, homoscedasticity and independence of the standardised innovations. Parsimonious parameterisation is based on log-likelihood ratio tests or on information criteria (e.g., AIC or BIC). The outcomes of such tests, however, depend on the particular point estimates of hyperparameters rather than on their entire distributions. Simulated distributions of the hyperparameter estimators, obtained with the Monte-Carlo simulation described in Section 4, give additional insight into the adequacy of the STS model. The simulated distributions give an indication as to whether or not the model tends to be overspecified in the sense that some state variables may be modelled as time invariant.

This study considers four models that differ in numbers of hyperparameters to be estimated with the ML method. The most complete model - Model 1 - is the one currently in use at Statistics Netherlands, but with the white noise component ${\epsilon }_t$ removed from the true population parameter ${\xi }_t.$ This component has turned out to have an implausibly large variance and disturbed estimation of other marginally significant hyperparameters (the seasonal and RGB disturbance variances) in the case of the DLFS. Removing the irregular component ${\epsilon }_t$ from the model has mitigated the instability in the two above-mentioned hyperparameters. This formulation implies that the population parameter ${\xi }_t$ does not exhibit irregularities that cannot be picked up by the stochastic structure of the trend and seasonal components. This assumption can be advocated by a relative rigidity of labour markets. Alterations of unemployment levels are usually gradual and therefore must be largely incorporated into the stochastic trend movements. The other three models are special cases of Model 1, i.e., all with the irregular component ${\epsilon }_t$ removed (see Table 5.1).

The simulated distributions of the hyperparameter estimators under Model 1 indicate that variance hyperparameters of the seasonal and, in particular, of the RGB component are often estimated to be close to zero. This causes bi-modality in the distribution of these variance estimates with a significant mass concentrated close to zero. Apart from that, an attempt to estimate both $\text{ln}\left({\widehat{\sigma }}_{\omega }^2\right)$ and $\text{ln}\left({\widehat{\sigma }}_{{\eta }_{\lambda }}^2\right),$ as in Model 1, distorts the other hyperparameters’ ML estimator distribution that is expected to be normal. For instance, normality in $\text{ln}\left({\widehat{\sigma }}_{v_3}^2\right),\text{ln}\left({\widehat{\sigma }}_{v_4}^2\right)$ and $\text{ln}\left({\widehat{\sigma }}_{v_5}^2\right)$ is severely violated with extreme outliers and/or a huge kurtosis (see Figure A.1 in Appendix, where the $\text{x}-$ axis is extended due to the outliers), while the corresponding variances are less likely to exhibit extreme values as they are supposed to fluctuate around 1. Making the seasonal component time-invariant, as in Model 2, hardly changes the situation for the trend and RGB hyperparameters. Instead, it may even be seen as suboptimal due to more extreme outliers and excess kurtosis in the distribution of all the five survey error hyperparameters (Figure A.2). By contrast, under both models where the $\text{RGB}-$ component is fixed over time (Models 3 and 4), all hyperparameter estimates corresponding to the survey error component have turned out to be normally distributed, see Figure A.3 and Figure A.4. Under Model 3, distributions are still skewed for the slope and seasonal components (skewness of -0.88 and -0.72, and kurtosis of 5.56 and 4.61, respectively). Fixing the seasonal hyperparameter to zero under Model 4 results in only a marginal improvement: the distribution of $\text{ln}\left({\widehat{\sigma }}_{{\eta }_R}^2\right)$ is negatively skewed (-0.81) with an excess kurtosis of 1.76.

This simulation evidence suggests that the preference in modelling the DLFS series may be given to the more parsimonious Model 3, where only the RGB disturbance variance is set equal to zero. However, since the RGB itself depends on the numbers of unemployed, its variance hyperparameter is retained for production purposes at Statistics Netherlands to secure sufficient flexibility against gradual changes in the underlying process.

The likelihood ratio test can be used to test if the hyperparameters of the seasonal and RGB components are significantly different form zero, since Models 2 through 4 are nested in Model 1. The test-statistic has very low values for all the three alternative models with respect to Model 1 (0, 0.18 and 0.18 for Models 2, 3 and 4, respectively, where the absence of differences between Models 2 and 1, as well as between Models 3 and 4 is due to a very low hyperparameter value of the seasonal component). These tests, thus, do not indicate that the more parsimonious models perform worse compared to Model 1. Another way of evaluating the adequacy of the four models is to compare their predictive power using the Root Mean Squared Differences (RMSD) between the GREG estimates and the one-step-ahead predictions for the signals. This is done for each wave separately: ${\text{RMSD}}^j\mathrm{<mo>=</mo>}1/\left(T-d\right){\displaystyle {\sum }_{t\mathrm{<mo>=</mo>}d}^T}{\left({\widehat{l}}_{t|t-1}^j-Y_t^j\right)}^2,$ with $d$ taken equal to 20, 30 and 60 months. Results presented in the Appendix (Table B.1), however, show that there is hardly any difference in the performance of the four models when applied to the original series. The more parsimonious models show a slight increase in the RMSD.

The distribution of the estimator of the survey error autoregressive parameter $\rho$ across the 1,000 simulated series does not seem to be affected by model reformulations: it approaches the normal distribution quite closely and ranges between 0 and 0.4 when $T\mathrm{<mo>=</mo>\; 114},$ which is in line with the approximation of its asymptotic distribution mentioned in Subsection 3.3. The range is slightly wider for the shorter time series and narrower when $T\mathrm{<mo>=</mo>\; 200.}$ The simulation procedure described in the previous section and the analysis of bootstrap methods that follows is performed separately for all the four models.

5.2 MSE estimation

The focus of this simulation study is MSE estimation for the trend and for the population signal, the latter being the sum of the trend and seasonal components. The performance of the Kalman filter and of the five MSE estimation methods discussed in Section 3 is evaluated by use of the relative bias and MSE of the MSE estimators. First, the filtered MSE estimates from (3.3), (3.4) and (3.7) are averaged over 1,000 simulations (where the average is denoted with a bar: ${\overline{\text{MSE}}}_{t|t})$ whereas the Kalman filter MSE estimates are averaged over 10,000 simulations, as mentioned at the beginning of Section 4. These averaged filtered MSE estimates for Model 3 (except for the $\text{AA}-$ method, see below why) are depicted in Figure 5.1 $–$ 5.4 for $T\mathrm{<mo>=</mo>\; 48},$ $T\mathrm{<mo>=</mo>\; 80},$ $T\mathrm{<mo>=</mo>\; 114},$ and $T\mathrm{<mo>=</mo>\; 200},$ respectively, skipping the first $d\mathrm{<mo>=</mo>\; 30}$ time points of the sample $(d$ should exceed the number of time points at the beginning of the series required to eliminate the effect of the diffuse filter initialization). Note that the analysis is based on filtered, rather than smoothed estimates, because filtered estimates better mimic the process of official figures production. MSEs in the four figures exhibit declining patterns, as expected, since the accuracy of the filtered estimates increases if more information becomes available over time for estimating the state variables. An exception is the true MSEs in Figure 5.2. A possible explanation is that, in this application, the signal MSEs are proportional to the signals themselves through the design-based standard errors, with the true MSEs being based on another (much larger) set of simulated series (50,000 for true MSEs; 1,000 for $\text{MSE}-$ estimators). Note that the lines in Figure 5.1 look much smoother because they are stretched over a smaller number of time points. Further, the patterns in Figure 5.2 $–$ 5.3 look more erratic because the scale of the $\text{y}-$ axis is finer, compared to Figure 5.1 and Figure 5.4.

The percentage relative bias is calculated as ${\text{RB}}_t^f\mathrm{<mo>=</mo>\; 100\; <mi>\%</mi>}\left({\overline{\text{MSE}}}_{t|t}^f/{\text{MSE}}_{t|t}^{\text{true}}-1\right),$ where $f$ defines a particular estimation method and ${\text{MSE}}_{t|t}^{\text{true}}$ is defined in (4.2). The percentage relative MSE biases averaged over time (skipping the first $d\mathrm{<mo>=</mo>\; 30}$ time points) for the signal, the trend and seasonal components are presented in Tables 5.2, 5.3, 5.4, and 5.5.

The main conclusions from the simulation study are as follows:

1. For $T\mathrm{<mo>=</mo>\; 48},$ and when averaged over time (starting from $t\mathrm{<mo>=</mo>\; 31}),$ the relative bias of the signal MSE obtained with the use of the Kalman filter is around -7 percent. This bias tends to decrease as the series length increases. The Kalman filter (KF) bias is quite small for the case of $T\mathrm{<mo>=</mo>\; 200},$ such that none of the estimation methods offers an improvement over the $\text{KF}-$ based $\text{MSE}-$ estimates. One could still apply the best estimation method with positive biases in order to get a range of values containing the true MSE.

2. The $\text{AA}-$ method turned out to be inapplicable to the models with marginally significant hyperparameters. When some of the hyperparameters are estimated close to zero, the matrix $I^{-1}\left({\widehat{\theta }}_{\sigma }^{\text{ML}}|{\rho }^a\right)$ is numerically either singular, leading to a failure in the procedure, or nearly singular. In the latter case, the asymptotic variance becomes excessively large and thus not reliable. Taking this into account, the $\text{AA}-$ method could only be considered for Model 4. As expected, the method performs poorly in short series, with positive biases of about 15 percent. The performance for $T\mathrm{<mo>=</mo>\; 114}$ and $T\mathrm{<mo>=</mo>\; 200}$ is comparable to that of the $\text{PT1}-$ bootstrap method, but significantly worse than the PT2 method’s performance.

3. As can be immediately observed, the use of the $\text{RR}-$ bootstrap results in a negative bias, whereas the use of the $\text{PT}-$ method produces a positive bias. Contrary to the claim of Rodriguez and Ruiz (2012) that their approach has better finite sample properties compared to the approach of Pfeffermann and Tiller (2005), the case of the DLFS suggests that the $\text{RR}-$ based MSE estimates, both the parametric and non-parametric ones, have larger negative biases than the $\text{KF}-$ based MSE estimates across all the models and series lengths (except for RR2 in Model 4 when $T\mathrm{<mo>=</mo>\; 48},$ and in Model 1 when $T\mathrm{<mo>=</mo>\; 80}$ and $T\mathrm{<mo>=</mo>\; 114}).$ While the $\text{PT}-$ bootstrap method is shown to have satisfactory asymptotic properties in Pfeffermann and Tiller (2005), Rodriguez and Ruiz (2012) illustrate the superiority of their method in small samples based on a simple model (a random walk plus noise). The present simulation study reveals that the $\text{RR}-$ method may not behave well in more complex applications. The $\text{PT}-$ methods have never produced negative biases for the DLFS, which makes these methods conservative (except for PT2 in Model 4 when $T\mathrm{<mo>=</mo>\; 48},$ with the negative bias still being smaller than that of the Kalman filter). Another striking outcome for $T\mathrm{<mo>=</mo>\; 48}$ is that the PT2 positive bias and the RR negative bias take on very large values in Model 3. However, with such a short series length and with so many non-stationary components like in the DLFS model, it is difficult to obtain reliable estimates from non-parametric bootstrap methods, since the burn-in period (or the diffuse sample) necessary for the non-parametric generation of the series takes more than a quarter of the series length (13 months out of 48).

4. For the series of lengths $T\mathrm{<mo>=</mo>\; 114}$ and $T\mathrm{<mo>=</mo>\; 80},$ the positive biases produced by the $\text{PT}2-$ method slightly exceed the $\text{KF}-$ biases in absolute value in models with insignificant hyperparameters (Models 1 and 2). In the more stable models (Models 3 and 4), the positive biases are smaller than the KF negative biases in absolute value. For $T\mathrm{<mo>=</mo>\; 48},$ bootstrap results are presented only for Models 3 and 4 (Models 1 and 2 that tend to be overspecified are not considered due to numerical problems). As could be expected, the biases are larger for this series length: the negative KF and RR biases become larger in absolute value, and so do the PT positive biases, with an exception of the above-mentioned result for PT2 in Model 4.

The signal MSE of Model 3, which could be considered as the better option for the production of official DLFS figures, is best estimated by the PT2 approach, with the relative bias of 1.4 and 1.9 percent for $T\mathrm{<mo>=</mo>\; 80}$ and $T\mathrm{<mo>=</mo>\; 114},$ respectively. The $\text{PT}2-$ bootstrap method also seems to be the best method for $T\mathrm{<mo>=</mo>\; 200},$ but, as already noted, the negative KF biases are already quite small for series of this length. For very short series like $T\mathrm{<mo>=</mo>\; 48},$ the parametric PT1 bootstrap seems to be the best option.

5. For both the $\text{PT}-$ and $\text{RR}-$ methods (except for RR2 in Model 4, $T\mathrm{<mo>=</mo>\; 48}),$ the absolute values of the relative biases are smaller in the case of the non-parametric approaches, compared to their parametric counterparts. The superiority of the non-parametric approach over the parametric one can be explained by the distorted normality of the error distribution in the models. Therefore, non-parametric bootstraps should be preferred unless time series are very short.

6. Apart from the bias of the MSE estimators, their variability may also give important insights into their reliability. To our knowledge, this has not been yet presented in the statistical literature. Tables 5.6 and 5.7 contain variances and MSEs of the four bootstrap MSE estimators for the signal, trend, and seasonal components for the two most interesting series lengths: $T\mathrm{<mo>=</mo>\; 48}$ and $T\mathrm{<mo>=</mo>\; 80}$ months (Models 1 and 2, as well as the asymptotic approximation, are not considered due to the aforementioned numerical problems). For both Model 3 and Model 4, the MSEs of the two PT MSE estimators are larger than the MSEs of the two RR MSE estimators. The RR MSE estimators’ seemingly superior performance, reflected by their smaller MSEs, is due to their smaller variances. The biases, however, are sometimes large enough to bring MSEs of these MSE estimators almost to the level of MSEs of the PT estimators. More importantly, the biases of the RR MSE estimators are mostly negative, often exceeding those of the Kalman filter. This phenomenon makes $\text{RR}-$ bootstraps hardly applicable in this application.

Apart from the above-mentioned simulation results, it is also interesting to see if the STS model-based approach still offers more precise predictors than the design-based variance estimates even after correcting for the hyperparameter uncertainty. For this purpose, $\text{STS}-$ model-based Root MSEs (RMSEs) obtained with the different MSE estimation procedures for the original series $\left(T\mathrm{<mo>=</mo>\; 114}\right)$ are compared to the standard errors (SEs) of the GREG estimator. Such Mean Differences in the Standard Errors (MDSE) under the time series model $m$ $\left(m\mathrm{<mo>=</mo>}\left\{\mathrm{1,2,3,4}\right\}\right)$ are defined as: ${\text{MDSE}}_m^f\mathrm{<mo>=</mo>}\mathrm{100\; <mi>\%</mi>}/\left(T-d\right){\displaystyle {\sum }_{t\mathrm{<mo>=</mo>}d}^T}\left[{\text{RMSE}}^f\left({\widehat{l}}_{t|t}^m\right)-\text{SE}\left(Y_t\right)\right]/\text{SE}\left(Y_t\right)$ and are presented in Table 5.8, with ${\widehat{l}}_{t|t}^m$ being the filtered estimate for the true population parameter, defined as trend plus seasonal, under model $m.$ Results are shown for the Kalman filter (labelled as “KF” in the table), i.e., when the hyperparameter uncertainty is neglected, as well as for cases when the five MSE estimation methods are applied to take the hyperparameter uncertainty into account. The true RMSEs from (4.2) are also compared to the GREG standard errors (see row “True” in Table 5.8). Note that the RGB and, particularly, the seasonal hyperparameter estimates obtained from the original DLFS data set are quite small. Therefore, there are no noticeable differences between the signal point-estimates of the four models. The AA, being the most unreliable approach, produces overestimated SEs (compare the 18- to 20-percent reduction based on the true RMSEs) due to nearly singular information matrices of the hyperparameter ML estimates. Keeping that in mind, one should feel more confident with the use of the $\text{PT}-$ estimators. Although the simulation study presented in this paper shows that PT2 usually outperforms the PT1 parametric approach, for this particular series, the $\text{PT1}-$ based SEs are closest to the true RMSEs, offering about a 20 percent reduction in the estimated GREG standard errors. This means that the model-based approach offers a significant variance reduction compared to the traditional design-based approach, even after accounting for the hyperparameter uncertainty.

• Previous
• Article main page
• Next