Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
112,486 characters · 8 sections · 70 citation commands
A dynamic factor model approach to incorporate Big Data in state space models for official statistics
There is an increasing interest among national statistical institutes (NSIs) to use data that are generated as a by-product of processes not directly related to statistical production purposes in the production of official statistics. Such data sources are sometimes referred to as “Big Data”; examples are time and location of network activity available from mobile phone companies, social media messages from Twitter and Facebook, sensor data, and internet search behaviour from Google Trends. A common problem with this type of data sources is that they are likely selective with respect to an intended target population. If such data sources are directly used to produce statistical information, then the potential selection bias of these data sources must be accounted for, which is often a hard task since Big Data sources are often noisy and generally contain no auxiliary variables, which are required for bias correction. These problems can be circumvented by using them as covariates in model-based inference procedures to make precise detailed and timely survey estimates, since they come at a high frequency and are therefore very timely. These techniques are known in the literature as small area estimation and nowcasting RaoMolina2015.
Official statistics are generally based on repeated samples. Therefore multivariate time series models are potentially fruitful to improve the precision and timeliness of domain estimates with survey data obtained in preceding reference periods and other domains. The predictive power of these models can be further improved by incorporating auxiliary series that are related with the target series observed with a repeated survey.
In this paper we investigate how auxiliary series derived from big data sources and registers can be combined with time series observed with repeated samples in high dimensional multivariate structural time series (STS) models. We consider Google Trends and claimant counts as auxiliary series for monthly unemployment estimates observed with a continuously conducted sample survey. Big Data sources have the problem that they are noisy and potentially (partly) irrelevant, and, as such, care must be taken when using them for the production of official statistics. We show that, by using a dynamic factor model in state space form, relevant information can be extracted from such auxiliary high-dimensional data sources, while guarding against the inclusion of irrelevant data.
Statistical information about a country's labour force is generally obtained from labour force surveys, since the required information is not available from registrations or other administrative data sources. The Dutch labour force survey (LFS) is based on a rotating panel design, where monthly household samples are observed five times with quarterly intervals. These figures are, however, considered too volatile to produce sufficiently reliable monthly estimates for the employed and the unemployed labour force at monthly frequency. For this reason Statistics Netherlands estimates monthly unemployment figures, together with its change, as unobserved components in a state space model where the observed series come from the monthly Dutch LFS, using a model originally proposed by Pfeffermann1991. This method improves the precision of the monthly estimates for unemployment with sample information from previous periods, and can therefore be seen as a form of small area estimation. In addition it accounts for rotation group bias bailar1975, serial correlation due to partial sample overlap, and discontinuities due to several major survey redesigns VanDenBrakel2015.
Time series estimates for the unemployment can be further improved by including related auxiliary series. The purpose is twofold. First, auxiliary series can further improve the precision of the time series predictions. In this regard, harvey2000 propose a bivariate state space model to combine a univariate series of the monthly unemployed labour force derived from the UK LFS, with the univariate auxiliary series of claimant counts. The latter series represents the number of people claiming unemployment benefits. It is an administrative source, which is not available for every country, and, as for the Netherlands, it can be affected by the same publication delay of the labour force series. Second, auxiliary series derived from Big Data sources like Google Trends are generally available at a higher frequency than the monthly series of the LFS. Combining both series in a time series model allows to make early predictions for the survey outcomes in real-time at the moment that the outcomes for the auxiliary series are available, but the survey data not yet, which is in the literature known as nowcasting, in other words, “forecasting the present".
In this paper, we extend the state space model used by Statistics Netherlands in order to combine the survey data with the claimant counts and the high-dimensional auxiliary series of Google Trends about job-search and economic uncertainty, as they could yield more information than a univariate one, which is not affected by publication lags and that can eventually be observed at a higher frequency than the labour force series.
This paper contributes to the existing literature by proposing a method to include a high-dimensional auxiliary series in a state space model in order to improve the (real-time) estimation of unobserved components. The model accounts for the rotating panel design underlying the sample survey series, combines series observed at different frequencies, and deals with missing observations at the end of the sample due to publication delays. It handles the curse of dimensionality that arises from including a large number of series related to the unobserved components, by extracting their common factors.
Besides claimant counts, the majority of the information related to unemployment is nowadays available on the internet; from job advertisements to resum\'e's templates and websites of recruitment agencies. We therefore follow the idea originating in Choi2009, askitas2009 and Suhoy2009 of using terms related to job and economic uncertainty, searched on Google in the Netherlands. Since 2004, these time series are freely downloadable in real-time from the Google Trends tool, on a monthly or higher frequency. As from the onset it is unclear which search terms are relevant, and if so, to which extent, care must be taken not to model spurious relationships with regards to the labour force series of interest, which could have a detrimental effect on the estimation of unemployment, such as happened for the widely publicized case of Google Flu Trends Lazer2014.
Our method allows to exploit the high-frequency and/or real-time information of the auxiliary series, and to use it in order to nowcast the unemployment, before the publication of labour force data. As the number of search terms related to unemployment can easily become large, we employ the two-step estimator of Doz2011, which combines factor models with the Kalman filter, to deal both with the high-dimensionality of the auxiliary series, and with the estimation of the state space model. The above-mentioned estimator is generally used to improve the nowcast of variables that are observed such as GDP (see giannone2008 and Hindrayanto2016 for applications to the US and the euro area), which is not the case for the unemployment. Nonetheless, DAmuriMarcucci2017, Naccaratoetal2018, and Maas2019 are all recent studies that use Google Trends to nowcast and forecast the unemployment, by treating the latter as known dependent variable in time series models where the Google searches are part of the explanatory variables. To the best of our knowledge, our paper is the first one to use Google Trends in order to nowcast the unemployment in a model setting that treats the latter variable as unobserved.
We evaluate the performance of our proposed method via Monte Carlo simulations and find that our method can yield large improvements in terms of Mean Squared Forecast Error $(\operatorname{MSFE})$ of the unobserved components' nowcasts. We then assess whether the accuracy of the unemployment's estimation and nowcast improves with our high-dimensional state space model, respectively from in-sample and out-of-sample results. The latter consists of a recursive nowcast. We do not venture into forecasting exercises as Google Trends are considered to be more helpful in predicting the present rather than the future of economic activities ChoiVarian2012. We conclude that Google Trends can significantly improve the fit of the model, although the magnitude of these improvements is sensitive to aspects of the data and the model specification, such as the frequency of observation of the Google Trends, the number of Google Trends' factors included in the model, and the level of estimation accuracy provided by the first step of the two-step estimation procedure.
The remainder of the paper is organized as follows. Section (ref) discusses the data used in the empirical analysis. Section (ref) describes the state space model that is currently used by Statistics Netherlands to estimate the unemployment. Section (ref) focuses on our proposed method to include a high-dimensional auxiliary series in the aforementioned model. Sections (ref) and (ref) report, respectively, the simulation and empirical results for our method. Section (ref) concludes.
The Dutch LFS is conducted as follows. Each month a stratified two-stage cluster design of addresses is selected. Strata are formed by geographical regions. Municipalities are considered as primary sampling units and addresses as secondary sampling units. All households residing on an address are included in the sample with a maximum of three (in the Netherlands there is generally one household per address). All household members with age of 16 or older are interviewed. Since October 1999, the LFS has been conducted as a rotating panel design. Each month a new sample, drawn according to the above-mentioned design, enters the panel and is interviewed five times at quarterly intervals. The sample that is interviewed for the $j^{th}$ time is called the $j^{th}$ wave of the panel, $j=1,\dots,5$. After the fifth interview, the sample of households leaves the panel. This rotation design implies that in each month five independent samples are observed. The generalized regression (GREG, i.e., design-based) estimator sarndal1992 is used to obtain five independent direct estimates for the unemployed labour force, which is defined as a population total. This generates over time a five-dimensional time series of the unemployed labour force. Table (ref) provides a visualization for the rotation panel design of the Dutch LFS.
Rotating panel designs generally suffer from Rotation Group Bias (RGB), which refers to the phenomena that there are systematic differences among the observations in the subsequent waves bailar1975. In the Dutch LFS the estimates for the unemployment based on the first wave are indeed systematically larger compared to the estimates based on the follow-up waves VanDenBrakel2015. This is the net results of different factors:
The Dutch labour force is subject to a one-month publication delay, which means that the sample estimates for month $t$ become available in month $t+1$. In order to have more timely and precise estimates of the unemployment, we extend the model by including, respectively, auxiliary series of weekly/monthly Google Trends about job-search and economic uncertainty, and monthly claimant counts, in the Netherlands.
Claimant counts are the number of registered people that receive unemployment benefits. The claimant counts for month $t$ become available in month $t+1$.
Google Trends are indexes of search activity. Each index measures the fraction of queries that include the term in question in the chosen geography at a particular time, relative to the total number of queries at that time. The maximum value of the index is set to be 100. According to the length of the selected period, the data can be downloaded at either monthly, weekly, or higher frequencies. The series are standardized according to the chosen period and their values can therefore vary according to the period's length Stephens-Davidowitz2015a. We use weekly and monthly Google Trends for each search term. Google Trends are available in real-time (i.e., they are available in period $t$ for period $t$, independently on whether the period is a week or a month).
The list of Google search terms used in the empirical analysis of this paper, together with their translation/explanation, is reported in Tables (ref) and (ref). A first set of terms (which is the one used in a previous version of this paper) was chosen by thinking of queries that could be made by unemployed people in the Netherlands. The rest of the terms has been chosen by using the Google Correlate tool and selecting the queries that are highly correlated to each term of the initial set, and that have a meaningful relation to unemployment and, more generally, economic uncertainty\footnote{Later in the paper we mention that we need non-stationary (e.g., persistent) Google Trends for our model. Correlations between non-stationary series can be spurious, and in this respect Google Correlate is not an ideal tool in order to choose search terms. In section (ref) we explain how to circumvent this problem.}.
Figure (ref) displays the time series of the five waves of the unemployed labour force, together with the claimant counts and an example of job-related Google query. They all seem to be following the same trend, which already shows the potential of using this auxiliary information in estimating the unemployment.
We first describe the model in use at Statistics Netherlands in Section 3.1. Next we explain how high-dimensional auxiliary series can be added to this model in Section 3.2.
The monthly sample size of the Dutch LFS is too small to produce sufficiently precise estimates directly. In the past, rolling quarterly figures were published on a monthly frequency. This has the obvious drawback that published figures are unnecessarily delayed since the reference period is the mid month of the rolling quarter. Also, real monthly seasonal effects are smoothed over the rolling quarter. Another problem that arose after the change from a cross-sectional survey to a rotating panel design in 2000, was that the effects of RGB became visible in the labour force figures. Both problems are solved with a structural time series (STS) model, that is used by Statistics Netherlands, since 2010, for the production of monthly statistics about the Dutch labour force VanDenBrakel2015. In a STS model, an observed series is decomposed in several unobserved components, such as a trend, a seasonal component, one or more cycles with a period longer than one year, regression components, and a white noise component. After writing an STS model in the state space form, the Kalman filter can be applied in order to estimate the unobserved components. See durbinkoopman2012 for an introduction to STS modelling.
Let $y^k_{j,t}$ denote the GREG estimate for the unemployment in month $t$ based on the sample observed in wave $j$. Now $\bm y^k_t = (y^k_{1,t}, \dots, y^k_{5,t})$ denotes the vector with the five GREG estimates for the unemployment in month $t$. The $y^k_{j,t}$ are treated as five distinct time series in a five dimensional time series model in order to account for the rotation group bias. The superscript $k > 1$ indicates that the vector is observed at the low frequency. We need this notation banbura2013 to distinguish between series observed at different frequencies, because later on we will make use of Google Trends which are available on a weekly basis. If $\bm y^k_t$ is observed at the monthly frequency, as in the case of the unemployed labour force, then $k=4,5$ if the high frequency series is observed at the weekly frequency, since a month can have either 4 or 5 weeks.
The unemployment is estimated, with the Kalman filter, as a state variable in a state space model where $\bm y^k_t$ represents the observed series. The measurement equation takes the form Pfeffermann1991, Brakel2009:
where $\bm\imath_5$ is a 5-dimensional vector of ones, and $\theta^{k,y}_t$, i.e. the unemployment, is the common population parameter among the five-dimensional waves of the unemployed labour force. It is composed of the level of a trend ($L_t$) and a seasonal component ($S_t$):
The transition equations for the level ($L_t$) and the slope ($R_t$) of the trend are, respectively:
which characterize a smooth trend model. This implies that the level of the trend is integrated of order 2, denoted as $I(2)$, which means that the series of the level is stationary (i.e., mean-reverting) after taking two times successive differences. The slope of the trend, $R^{k,y}_t$, is a first-order integrated series, denoted as $I(1)$. This state variable represents the change in the level of the trend, $L^{k,y}_t$, and not in the unemployment, $\theta^{k,y}_t$, directly. Nevertheless, since the $I(2)$ property of the unemployment is driven by its trend, and not by its seasonal component, the change in $\theta^{k,y}_t$ will also mainly be captured by $R^{k,y}_t$, and we can therefore consider the latter as a proxy for the change in unemployment. The model originally contained an innovation term for the population parameter $\theta^{k,y}_t$. However, the maximum likelihood estimate for its variance tended to be zero and bollineni2017 showed via simulations that it is better to not include this term in the model.
The trigonometric stochastic seasonal component allows for the seasonality to vary over time, and it is modeled as in durbinkoopman2012:
where $h_l = \frac{\pi l}{6}$, for $l=1,\dots,6$.
The second component in equation (ref), $ \bm \lambda^k_t = (\lambda^k_{1,t}, \ldots, \lambda^k_{5,t})^t$, accounts for the RGB. Based on the factors that contribute to the RGB, as mentioned in Section (ref), the response observed in the first wave is assumed to be the most reliable one and not to be affected by the RGB Brakel2009. Therefore it is assumed that $\lambda^k_{1,t} = 0$. The remaining four components in $\bm \lambda^k_t$ are random walks that capture time-dependent differences between the follow-up waves with respect to the first wave:
As a result the Kalman filter estimates for $\theta^{k,y}_t $ in (ref) are benchmarked to the level of the GREG series of the first wave.
The third component in equation (ref), $ \bm e^k_t = (e^k_{1,t}, \ldots, e^k_{5,t})^t$, models the autocorrelation among the survey errors ($e^k_{j,t}$) in the follow-up waves due to the sample overlap of the rotating panel design. In order to account for this autocorrelation, the survey errors are treated as state variables, which follow the transition equation below.
with $\widehat{ \operatorname{var}} \left(y^k_{j,t}\right)$ being the design variance of the GREG estimates $y^k_{j,t}$. The scaled sampling errors, $\tilde{e}^k_{j,t}$, for $j=1,\dots,5$, account for the serial autocorrelation induced by the sampling overlap of the rotating panel. Samples in the first wave are observed for the first time and therefore its survey errors are not autocorrelated with survey errors of previous periods. The survey errors of the second to fifth wave are correlated with the survey errors of the previous wave three months before. Based on the approach proposed by Pfeffermann1998, Brakel2009 motivate that these survey errors should be modelled as an AR(3) process, without including the first two lags. Moreover, the survey errors of all waves are assumed to be proportional to the standard error of the GREG estimates. In this way the model accounts for heterogeneity in the variances of the survey errors, which are caused by changing sample sizes over time. As a result the maximum likelihood estimates of the variances of the scaled sampling errors, $\sigma^2_{\nu_{j}}$, will have values approximately equal to one.
The structural time series model (ref) as well as the models proposed in the following sections are fitted with the Kalman filter after putting the model in state space form. We use an exact initialization for the initial values of the state variables of the sampling error, and a diffuse initialization for the other state variables. It is common to call hyperparameters the parameters that define the stochastic properties of the measurement equation and the transition equation of the state space model. These are the parameters that are assumed to be known in the Kalman filter durbinkoopman2012. In our case the hyperparameters are $\delta$ and all the parameters that enter the covariance matrices of the innovations. These hyperparameters are estimated by maximum likelihood using the Broyden-Fletcher-Goldfarh-Shanno (BFGS) optimization algorithm. The additional uncertainty of using maximum likelihood estimates for the hyperparameters in the Kalman filter is ignored in the standard errors of the filtered state variables. Since the observed time series contains 185 monthly periods, this additional uncertainty can be ignored. See also bollineni2017 for details. Both the simulation and estimation results in Sections (ref) and (ref) are obtained using the statistical software R.
Assuming normality of the innovations is common in state space models because the hyperparameters of the model are estimated by maximizing a Gaussian log-likelihood which is evaluated by the Kalman filter. Moreover, under normality, the Kalman filter yields the minimum variance unbiased estimator of the state variables. Nonetheless, as long as the state space model is linear, if the true distribution of the error terms is non-Gaussian, then the Kalman filter still provides the minimum variance linear unbiased estimator of the state variables durbinkoopman2012. In this case we can further rely on quasi maximum likelihood (QML) theory in order to perform inference based on the QML estimates of the hyperparameters. This means that the hyperparameters can still be consistently estimated by maximizing the Gaussian log-likelihood (or in general, as Gourierouxetal1984 argue, a density function that belongs to the family of linear exponential distributions), but we shall use, if needed, the appropriate expression for the covariance matrix of the QML estimators, which should capture the additional uncertainty caused by the model's misspecification hamilton1994. In Appendix (ref) we conduct a Monte Carlo simulations study and find that deviations from normality are not of concern for the performance our method.
This time series model addresses and solves the mentioned problems with small sample sizes and rotation group bias. Every month a filtered estimate for the trend ($ L^{k,y}_{t}$) and the population parameter, which is defined as the filtered trend plus the filtered seasonal effect ($ \theta^{k,y}_{t} = L^{k,y}_{t} + S^{k,y}_{t}$), are published in month $t+1$. The time series model uses sample information from previous months in order to obtain more stable estimates. The estimates account for RGB by benchmarking the estimates for $ L^{k,y}_{t}$ and $ \theta^{k,y}_{t}$ to the level of the first wave, which makes them comparable with the outcomes obtained under the cross-sectional design before 2000.
We now introduce some further notation to distinguish between in-sample estimates and out-of-sample forecasts. In the case of in-sample estimates, $\hat \theta^{k,y}_{t|\Omega_t}$ denotes the filtered estimate of the population parameter $\theta^{k,y}_{t}$, assuming that all data for time $t$ is released and available at time $t$. We therefore condition on the information set $\Omega_t$ which does not contain any missing data at time $t$. In the case of out-of-sample forecasts, we condition on the data set $\Omega_t^{-}$ that is actually available in real time at time $t$. For instance, $\bm y_t^k$ only gets published during moth $t+1$, and is therefore not available yet at time $t$, and not part of $\Omega_t^{-}$. Thus $\hat \theta^{k,y}_{t|\Omega_t^{-}}$ is the filtered forecast for $\theta^{k,y}_{t}$, based on the information that is available at time $t$. Under model ((ref)), which does not contain auxiliary information other than the labour force series, $\hat \theta^{k,y}_{t|\Omega_t^{-}}$ is in fact the one-step-ahead prediction $\hat \theta^{k,y}_{t|\Omega_{t-1}}$, since $\bm y_t^k$ is not available yet in month $t$, but $\bm y_{t-1}^k$ is; therefore, $\Omega_t^{-} = \Omega_{t-1} = \{ \bm y^k_{t-1}, \bm y^k_{t-2}, \ldots \}$.
To improve precision and timeliness of the monthly unemployment figures, we extend the labour force model by including auxiliary series of weekly/monthly Google Trends about job-search and economic uncertainty, and monthly claimant counts, in the Netherlands. Since the claimant counts for month $t$ become available in month $t+1$, it is anticipated that this auxiliary series is particularly useful to further improve the precision of the trend and population parameter estimates after finalizing the data collection for reference month $t$. The Google Trends come at a higher frequency already during the reference month $t$. It is therefore anticipated that these auxiliary series can be used to make first provisional estimates for the trend and the population parameter of the LFS during month $t$, when the sample estimates $\bm y^k_t$ are not available, but the Google Trends become available on weekly basis.
Weekly and monthly Google Trends are throughout the paper denoted by $\bm x^{GT}_t$ and $\bm x^{k,GT}_t$, respectively. We denote the dimension of the vector $\bm x^{GT}_t$ by $n$, which can be large. In addition, we can expect the Google Trends to be very noisy, such that the signal about unemployment contained in them is weak. We therefore need to address the high-dimensionality of these auxiliary series, in order to make the dimension of our state space model manageable for estimation, and extract the relevant information from these series. For this purpose we employ a factor model which achieves both by retaining the information of these time series in a few common factors.
Moreover, when dealing with mixed frequency variables and with publication delays, we can encounter “jagged edge" datasets, which have missing values at the end of the sample period. The Kalman filter computes a prediction for the unobserved components in presence of missing observations for the respective observable variables.
The two-step estimator by Doz2011 combines factor models with the Kalman filter and hence addresses both of these issues. In the remainder of this section we explain how this estimator can be employed to nowcast the lower-frequency unobserved components of the labour force model using information from higher-frequency or real-time auxiliary series.
We consider the following state space representation of the dynamic factor model for the Google Trends, with respective measurement and transition equations, as we would like to link it to the state space model used to estimate the unemployment (ref):
where $\bm x^{GT}_t$ is a $n\times1$ vector of observed series, $\bm f_t$ is a $r \times 1$ vector of latent factors with $r \ll n$, $\bm \varLambda$ is a $n \times r$ matrix of factor loadings, $\bm \varepsilon_t$ is the $n \times 1$ vector of idiosyncratic components and $\bm \varPsi$ its $n \times n$ covariance matrix; $\bm u_t$ is the $r\times1$ vector of factors' innovations and $\bm I_r$ is a $r\times r$ identity matrix (which follows from the identification conditions used in principal component analysis since the factors are only identified up to rotation). Notice that the dynamic equation for $\bm f_t$ implies that we are making the assumption that $\bm x^{GT}_t$ is $I(1)$ of dimension $n$, and $\bm f_t$ is $I(1)$ of dimension $r$. Later in this section the need of this assumption will become clearer; the intuition behind it is that the factors and the change in unemployment, $R^{k,y}_t$, must be of the same order of integration.
Among others, Bai2004 proves the consistency of the estimator of $I(1)$ factors by principal component analysis (PCA), under the assumptions of limited time and cross-sectional dependence and stationarity of the idiosyncratic components, $\bm \varepsilon_t$, and non-trivial contributions of the factors to the variance of $\bm x_t$.\footnote{For the exact formulation we refer to Assumptions A-D in Bai2004.} We assume no cointegrating relationships among the factors. We further assume normality of the innovations for the same reasons outlined in Section (ref).
The consistency of the two-step estimator has been originally proven in the stationary framework by Doz2011, and extended to the nonstationary case by Barigozzi2017.
In the first step, the factors ($\bm f_t$), the factor loadings ($\bm \varLambda$), and the covariance matrix of the idiosyncratic components ($\bm \varPsi$) in model (ref) are estimated by PCA as in Bai2004. The matrices $\bm \varLambda$ and $\bm \varPsi$ are then replaced, in model (ref), by their estimates $\hat{\bm \varLambda}$ and $\hat{\bm \varPsi}=\operatorname{diag}\left(\hat{\psi}_{11}, \dots , \hat{\psi}_{nn} \right)$ obtained in this first step. These estimates are kept as fixed in the second step, because their high-dimensionality and associated curse of dimensionality complicates re-estimation by maximum likelihood. Moreover, restricting the covariance matrix of the idiosyncratic components $\bm \varPsi$ as being diagonal is standard in the literature\footnote{The specification of the dynamic factor factor model with spherical idiosyncratic components is often called “approximate" dynamic factor model. Doz2011 and Barigozzi2017 mention that misspecifications of this model arising from time or cross-sectional dependence of the idiosyncratic components, do not affect the consistency of the two-step estimator of the unobserved common factors, if $n$ is large.}
In order to make use of the auxiliary series to nowcast the unemployment, we stack together the measurement equations for $\bm y^k_t$ and $\bm x^{k,GT}_t$, respectively (ref) and the first equation of (ref) with $\bm \varLambda$ and $\bm \varPsi$ replaced, respectively, by $\hat{\bm \varLambda}$ and $\hat{\bm \varPsi}$, and express them at the lowest frequency (in our case the monthly observation's frequency of $\bm y^k_t$). The transition equations for the RGB and survey error component in combination with the rotation scheme applied in the Dutch LFS hamper a formulation of the model on the high frequency. This means that $\bm x^{GT}_t$ needs to be first temporally aggregated from the high to the low frequency (either before or after the first step which estimates $\bm \varLambda$ and $\bm \varPsi$). Since $\bm x^{GT}_t$ are the $I(1)$ weekly Google Trends, which are flow variables as they measure the number of queries made during each week, they are aggregated according to the following rule banbura2013:
The aggregated $\bm x^{k,GT}_{j,t}$ are then rescaled in order to be bounded again between 0 and 100. The subscript $j$ allows for real-time updating of the aggregated Google Trends in week $j$ when new data become available. As such, this index indicates that we aggregate weeks 1 up to $j$. When $j=k$ we are at the end of the month, and we simply write $\bm x^{k,GT}_{t}$ to indicate the end-of-month aggregate value.
In order to get the final model, we also include a measurement equation for the univariate auxiliary series of the claimant counts, assuming that its state vector, $\theta^{k,CC}_t$, has the same composition of our population parameter $\theta^{k,y}_t$ (i.e., composed of a smooth trend and a seasonal component):
The last equality allows the innovations of the trends' slopes, $R^{k,y}_t$ and $R^{k,CC}_t$, and of the factors of the Google Trends, to be correlated. harvey2000 show that there can be potential gains in precision, in terms of Mean Squared Error $(\operatorname{MSE})$ of the Kalman filter estimators of $\theta^{k,y}_t$, $L^{k,y}_t$, and $R^{k,y}_t$, if the correlation parameters $|\rho|$s are large. Specifically, if $|\rho_{CC}| = 1$, then $\bm y^k_t$ and $x^{k,CC}_t$ have a common slope. This means that $\bm y^k_t$ and $x^{k,CC}_t$ are both $I(2)$, but there is a linear combination of their first differences which is stationary. Likewise, if $|\rho_{m,GT}| = 1$ then the $m^{\text{th}}$ factor of the Google Trends and the change in unemployment, $R^{k,y}_t$, are cointegrated (i.e., they have the same source of error). This is why we need the elements of the vector in (ref) to have the same order of integration, and it is via this correlation parameters that we exploit the auxiliary information.
The second step of the estimation procedure consists of estimating the remaining hyperparameters of the whole state space model (equations (ref)-(ref)) by maximum likelihood, and applying the Kalman filter to re-estimate $\bm f^k_t$ and to nowcast the variables of interest, $\theta^{k,y}_t$, $L^{k,y}_t$, and $R^{k,y}_t$, providing unemployment estimates in real-time before LFS data become available: $\hat{\theta}^{k,y}_{t|\Omega_t^{-}}$, $\hat{L}^{k,y}_{t|\Omega_t^{-}}$, and $\hat{R}^{k,y}_{t|\Omega_t^{-}}$ are the filtered nowcasts of, respectively, $\theta^{k,y}_t$, $L^{k,y}_t$, and $R^{k,y}_t$ based on the information set $\Omega_t^{-}$ available in month $t$. The information set in this case is $\Omega_t^{-} = \{\bm x^{k,GT}_{t}, \bm y^k_{t-1}, x^{k,CC}_{t-1}, \bm x^{k,GT}_{t-1}, \ldots \}$. Note that, contrary to Section (ref), we now talk about “nowcast" instead of “forecast" of $\theta^{k,y}_t$ because a part of the data (the Google Trends) used in model (ref)-(ref) is now available in month $t$.
Some remarks are in place. First, although in Section (ref) we mentioned that Statistics Netherlands publishes only $\hat{L}^{k,y}_t$ and $\hat{\theta}^{k,y}_t$ as official statistics for the unemployment, we are also interested in the estimation/nowcast accuracy of $R^{k,y}_t$ since it is the state variable of the labour force model that is directly related to the auxiliary series.
Second, note that in model (ref) we do not make use of the superscript $k$, meaning that the first step of the estimation can be performed on the high frequency (weekly in our empirical case) variables. Since in each week we can aggregate the weekly Google Trends to the monthly frequency, we can use the information available throughout the month to update the estimates of $\bm \varLambda$ and $\bm \varPsi$. If the correlations between the factors and the trend's slope of the target variable are large, this update should provide a more precise nowcast of $R^{k,y}_t$, $L^{k,y}_t$ and $\theta^{k,y}_t$.
Third, we allow the factors of the Google Trends to be correlated with the change in unemployment and not with its level for two reasons: first, a smooth trend model is assumed for the population parameter, which means that the level of its trend does not have an innovation term. Second, it is reasonable to assume that people start looking for a job on the internet when they become unemployed, and hence their search behaviour should reflect the change in unemployment rather than its level.
Fourth, while our method to include auxiliary information in a state space model is based on the approach proposed by harvey2000, the factors of the high-dimensional auxiliary series could also be included as regressors in the observation equation for the labour force. However, in such a model, the main part of the trend, $L_t^{k,y}$, will be explained by the auxiliary series in the regression component. As a result, the filtered estimates for $L_t^{k,y}$ will contain a residual trend instead of the trend of the unemployment. Since the filtered trend estimates are the most important target variables in the official monthly publications of the labour force, this approach is not further investigated in this paper.
Finally, we refer the reader to Appendices (ref), (ref), and (ref) for a detailed state space representation of the labour force model when, respectively, a univariate, a high-dimensional or both type of auxiliary series are included. We further refer to Appendices (ref) and (ref) for an illustration on how to include the lags of the factors and how to model their cycle or seasonality, within our proposed high-dimensional state space model.
We next conduct a Monte Carlo simulations study in order to elucidate to which extent our proposed method can provide gains in the nowcast accuracy of the unobserved components of interest. For this purpose, we consider a simpler model than the one used for the labour force survey. Here $y^k_t$ is univariate following a smooth trend model, and $\bm x^k_t$ represents the $(100 \times 1)$-dimensional auxiliary series with one common factor ($r = 1$).
We allow the slope and factor's innovations to be correlated, and we investigate the performance of the method for increasing values of the correlation parameter $\rho \in [0,0.2,0.4,0.6,0.8,0.9,0.99]$. The auxiliary variable $\bm x^k_t$ has the same frequency of $y^k_t$ and it is assumed that all $\bm x^k_t$ are released at the same time without publication delays. The nowcast is done concurrently, i.e. in real-time based on a recursive scheme. This means that in each time point of the out-of-sample period, the hyperparameters of the model are re-estimated by maximum likelihood, extending the same used up to that period. This is done in the third part of the sample, always assuming that $y^k_t$ is not available at time $t$, contrary to $\bm x^k_t$. This implies that the available data set in period $t$ equals $\Omega_t^{-} = \{\bm x^k_t, y^k_{t-1}, \bm x^k_{t-1}, y^k_{t-2}, \ldots \}$. The sample size is $T = 150$ and the number of simulations is $n_{\text{sim}} = 500$.
We consider three specifications for the idiosyncratic components and the factor loadings:
Let $\bm \alpha^k_t = \left(L^k_t, R^k_t, f^k_t \right)'$ denote the vector of state variables and $\hat \bm \alpha^k_{t|\Omega_t^{-}}$ its estimates based on the information available at time $t$. The results from the Monte Carlo simulations are shown in Table (ref). We always report the $\operatorname{MSFE}$, together with its variance and bias components, of the Kalman filter estimator of $\bm \alpha^k_t$, relative to the same measures calculated from the model that does not include the auxiliary series $\bm x^k_t$. Recall that the latter comes down to making one-step-ahead predictions.
where $h$ is the size of the of out-of-sample period.
In every setting, both the bias and the variance of the $\operatorname{MSFE}$ tend to decrease with the magnitude of the correlation parameter. The improvement is more pronounced for the slope rather than the level of the trend. For the largest value of the correlation, with respect to the model which does not include auxiliary information, the gain in $\operatorname{MSFE}$ for the level and the slope is, respectively, of around 25% and 75%. Moreover, for low values of $\rho$, the $\operatorname{MSFE}$ does not deteriorate with respect to the benchmark model. This implies that our proposed method is robust to the inclusion of auxiliary information that does not have predictive power for the state variables of interest. In Appendix (ref) we report and examine additional simulation results with non-Gaussian idiosyncratic components, and draw the same conclusions discussed above for the $\operatorname{MSFE}$ and the variance of the state variables' nowcasts. The bias instead worsens while deviating from Gaussianity, but it does not affect the $\operatorname{MSFE}$ as it only accounts for a small part of the latter measure. We therefore conclude that the performance of our method is overall robust to deviations from Gaussianity of the idiosyncratic components.
The decision to focus the simulation study on the nowcast (rather than the in-sample) performance of our method, is motivated by the fact that the added value of the Google Trends over the claimant counts is their real-time availability, which can be used to nowcast the unemployment. Nonetheless, for completeness, in the empirical application of the next section we report the results also for the in-sample performance of our method.
In this section we present and discuss the results of the empirical application of our method to nowcasting the Dutch unemployment using the auxiliary series of claimant counts and Google Trends related to job-search and economic uncertainty.
As explained in Section (ref), the Google series used in the model must be $I(1)$. We therefore test for nonstationarity in the Google Trends with the Elliott1996 augmented Dickey-Fuller (ADF) test, including a constant and a linear trend. We control for the false discovery rate as in Moon2012, who employ a moving block bootstrap approach that accounts for time and cross-sectional dependence among the units in the panel.
Before proceeding with the estimation of the model by only including the Google Trends that resulted as being $I(1)$ from the multiple hypotheses testing, we carry out an additional selection of the $I(1)$ Google Trends by “targeting" them as explained and motivated below.
BAING2008target point out that having more data to extract factors from is not always better. In particular, if series are added that have loadings of zero and are thus not influenced by the factors, these will make the estimation of factors and loadings by PCA deteriorate, as PCA assigns a non-zero weight to each series in calculating the estimated factor as a weighted average. BAING2008target recommend a simple strategy to filter out irrelevant series (in our case Google search terms) and improve the estimation of the factors, which they call “targeting the predictors". In this case an initial regression of the series of interest is performed on the high-dimensional input series to determine which series are (ir)relevant. The series that are found to be irrelevant are discarded and only the ones that are found to be relevant are kept to estimate the factors and loadings from. In particular, they recommend the use of the elastic net hastie2005, which is a penalized regression technique that performs estimation and variable selection at the same time by setting the coefficients of the irrelevant variables to 0 exactly. After performing the elastic net estimation, only the variables with non-zero coefficients are then kept. As we do not observe our series of interest directly, we need to adapt their procedure to our setting. To do so we approximate the unobserved unemployment by its estimation from the labour force model without auxiliary series. Specifically, we regress the differenced estimated change in unemployment from the labour force model without auxiliary series, $\Delta \hat{R}^{k,y}_t$, on the differenced $I(1)$ Google Trends using the elastic net penalized regression method, which solves the following minimization problem:
where
The tuning parameters $\lambda$ and $\alpha$ are selected from a two-dimensional grid in order to minimize the Schwarz1978 Bayesian information criterion (BIC). Notice that performing the penalized regression on the differenced (and therefore stationary) data, also allows us to avoid the inclusion in the model of Google Trends that have spurious relations with the change in unemployment.
We both consider estimating the final model with all Google Trends included and with only the selected Google Trends included, thereby allowing us to assess the empirical effects of targeting. The final number of nonstationary Google Trends included in the model, $n$, may differ depending on whether we use the weekly Google Trends aggregated to the monthly frequency according to equation (ref), or the monthly Google Trends. Whenever we apply PCA, the Google Trends are first differenced and standardized.
We further need to make sure that the stationarity assumption of the idiosyncratic components is maintained. Therefore, after having estimated the factors by PCA in model (ref), we test which of the idiosyncratic components $\bm \varepsilon_t$ are $I(1)$ with an ADF test without deterministic components, by controlling for multiple hypotheses testing as in Moon2012. The $I(1)$ idiosyncratic components are modelled as state variables in (ref), with the following transition equation:
with usual normality assumptions on the $\bm \xi^k_t$. The covariance matrix of the idiosyncratic components $\bm \varPsi$ is therefore estimated on the levels of the $I(0)$ idiosyncratic components and the first differences of the $I(1)$ idiosyncratic components. Appendix (ref) provides a toy example that elucidates the estimation procedure.
Finally, we notice that although the first step of the two-step estimation procedure is meant to avoid estimating $\bm \varPsi$ and $\bm \varLambda$ by maximum likelihood (since they are large matrices), this pre-estimation may affect the explanatory power of the Google Trends. We here propose two different ways to obtain (possibly) more accurate estimates of these two matrices:
Later in this section we check how sensitive our empirical results are to the different estimates of $\bm \varPsi$ and $\bm \varLambda$. For the second type of estimation method discussed above, we only perform one additional iteration of the two-step procedure due to its computational burden.
We present empirical results for the in-sample estimates and out-of-sample forecasts. With the in-sample estimates we evaluate to which extent the auxiliary series improve the precision of the published monthly unemployment estimates after finalizing the data collection. With the out-of-sample forecasts we evaluate to which extent the auxiliary series improve the precision of provisional estimates in a nowcast procedure during the period of data collection. We always estimate four different models: the labour force model without auxiliary series (baseline), the labour force model with auxiliary series of claimant counts (CC), of Google Trends (GT) and of both (CC & GT). We compare the latter three models to the baseline one with the in-sample and out-of-sample exercises. The period considered for the estimation starts in January 2004 and ends in May 2019 ($T=185$ months). The out-of-sample nowcasts are conducted in real-time (concurrently) in the last three years of the sample based on a recursive scheme: each week or month, depending on whether we use weekly or monthly Google Trends, the model, including its hyperparameters, is re-estimated on the enlarged sample now extended by the latest observations, while assuming that the current observations for the unemployed labour force and the claimant counts are missing. Analogously, when the Google Trends are first targeted with the elastic net, the targeting is re-executed in each week or month of the out-of-sample period on the updated sample.
We define the measure of in-sample estimation accuracy $\widehat{\operatorname{MSE}}(\hat{\bm \alpha}^k_{t|\Omega_t}) = \frac{1}{T-d} \sum_{t=d+1}^{T} \hat{\bm P}^k_{t|\Omega_t}$, where $\hat \bm \alpha^k_{t|\Omega_t}$ is the vector of Kalman filter estimates of the state variables, $\hat{\bm P}^k_{t|\Omega_t}$ is its estimated covariance matrix in month $t$, and $d$ is the number of state variables that are needed to estimated the labour force model without auxiliary series, and that need a diffuse initialization for their estimation ($d=17$). The measure of nowcast accuracy, $\widehat{\operatorname{MSFE}}(\hat{\bm \alpha}^k_{t|\Omega_t^{-}}) = \frac{1}{h} \sum_{t=T-h+1}^{T} \hat{\bm P}^k_{t|\Omega_t^{-}}$, is the average of the nowcasted covariance matrices in the $h$ prediction months. When weekly Google Trends are used, $ \hat{\bm P}^k_{t|\Omega_t^{-}} = \frac{1}{k} \sum_{j=1}^{k} \hat{\bm P}^k_{j|\Omega_{j,t}^{-}}$, where $\hat{\bm P}^k_{j|\Omega_{j,t}^{-}}$ is the nowcasted covariance matrix for the prediction in week $j$ of month $t$, and $\Omega_{j,t}^{-} = \{\bm x^{k,GT}_{j,t}, \bm y^k_{t-1}, x^{k,CC}_{t-1}, \bm x^{k,GT}_{t-1}, \ldots\}$ is in this case the available information set in week $j$ of month $t$. This is because the nowcast is done recursively throughout the weeks of the out-of-sample period. We always report the relative $\widehat{\text{MS(F)E}}$ with respect to the baseline model; values lower than one are in favour of our method. We note that nowcasting under the baseline model without auxiliary series and the baseline model extended with claimant counts comes down to making one-step-ahead predictions. Expressions for $\hat \bm \alpha^k_{t|\Omega_t}$, $\hat \bm \alpha^k_{t|\Omega_t^{-}}$ and their covariance matrices, $\hat{\bm P}^k_{t|\Omega_t}$ and $\hat{\bm P}^k_{t|\Omega_t^{-}}$, are given by the standard Kalman filter recursions, see e.g. durbinkoopman2012.
The initial values of the hyperparameters for the maximum likelihood estimation are equal to the estimates for the labour force model obtained in VanDenBrakel2015. We use a diffuse initialisation of the Kalman filter for all the state variables except for the 13 state variables that define the autocorrelation structure of the survey errors, for which we use the exact initialisation of bollineni2017.
We use the three panel information criteria proposed by baing2002 which we indicate, as in baing2002, with $IC_1$, $IC_2$ and $IC_3$, in order to choose how many factors of the Google Trends to include in the model\footnote{In this paper, if for instance the information criterion $IC_1$ suggests to include 2 factors, we indicate as $IC_1=2$.}. When the Google Trends are targeted with the elastic net, the information criteria suggest to include one or two factors. In the empirical analysis we check the sensitivity of the results with respect to these two different numbers of factors included in the model.
We employ a wilks1938 likelihood ratio (LR) test to assess whether the correlation parameters are significantly different from zero, and hence adding the auxiliary information might yield a significant improvement from the baseline model. Specifically, we indicate with $\rho_{CC} = 0$, $\rho_{1,GT} = 0$ and $\rho_{2,GT} = 0$ the null hypotheses for the individual insignificance of the correlation parameter with, respectively, the claimant counts, and the first and second factor (when present) of the Google Trends. With $\bm \rho_{GT} = \boldsymbol{0}$ and $\bm \rho = \boldsymbol{0}$ we instead indicate the null hypotheses for the joint insignificance of, respectively, the correlations with the Google Trends' factors, and all correlation parameters. If the true distribution of the error terms is non-Gaussian, the LR test, based on the QML estimates, does not generally keep having, under the null hypothesis, an asymptotic $\chi^2$ distribution with degrees of freedom equal to the number of restrictions. One exception is when the covariance matrix of the error terms from a regression involving observed variables, is replaced by a consistent estimator prior to the maximization of the log-likelihood GourierouxMonfort1993. In our case, if the idiosyncratic components of the Google Trends, $\bm \varepsilon^{k,GT}_t$, are the only error terms not being normally-distributed, we may fall into this exception. The covariance matrix $\bm \varPsi$ is indeed replaced, for the maximization of the log-likelihood, by its consistent PCA estimator obtained in the first step of the two-step estimation procedure. Nonetheless, in the setting of GourierouxMonfort1993 the regressors are observed, whereas in our case the latter are the unobserved factors. Consequently, it is not trivial to asses whether our model specification indeed falls into the above-mentioned exception. A formal proof for this is beyond the scope of this paper, but in Appendix (ref) we conduct a simulation study in order to obtain the finite-sample probability density of the LR test under misspecifications of the distribution of the idiosyncratic components. We conclude that the distribution of the LR test is not affected by these misspecifications. At the end of this section we show that that there is no evidence that the error terms other than $\bm \varepsilon^{k,GT}_t$, are not normally-distributed. We should therefore be able perform inference based on the usual asymptotic distribution of the LR test.
Table (ref) reports the estimated hyperparameters for the four models, as well as the respective value for the maximized log-likelihood, the relative measures of in and out-of-sample performance, and the p-values from the LR tests, when the monthly Google Trends are used.
The maximum likelihood estimates for the standard error of the seasonal components' disturbance terms tend to zero, indicating that the seasonal effects are time invariant.
Recall from equation ((ref)) that the variances of the scaled sampling errors, $\sigma^2_{\nu_{j}}$, should take values close to one. Their estimates are divided by (1 - $\hat{\delta}^2$) and are always slightly larger than one, which is an indication that the variance estimates of the GREG estimates, used to scale the sampling errors in equation ((ref)), somewhat underestimate the real variance of the GREG estimates.
The correlation with the claimant counts is estimated to be above 0.9, and remains large and significant when including the Google Trends. Similar conclusions can be drawn for the correlations with the Google Trends' factors, when the Google Trends are targeted with the elastic net, and 39 of them are included in the model. When the additional targeting is not applied, and the 162 $I(1)$ Google Trends are directly included in the model, the correlation parameter with the first factor of the Google Trends is instead always small and insignificant (in this setting we do not include more than one factor). Moreover, for the same number of factors, targeting the Google Trends always yield a better performance in terms of estimation and nowcast accuracy of the state variables of interest, with respect to not targeting them. For this reason, we focus the remaining analysis of the empirical results only on the targeted Google Trends.
The best results in terms of both estimation and nowcast accuracy of all the state variables, is achieved by the CC & GT model with one factor, yielding a gain of, respectively, around 40% and 20% for $\hat{R}^{k,y}_{t|\Omega_t^{-}}$, and around 20% and 25% for both $\hat{L}^{k,y}_{t|\Omega_t^{-}}$ and $\hat{\theta}^{k,y}_{t|\Omega_t^{-}}$, with respect to the baseline model. Note that this implies that the above-mentioned model outperforms also the model that contains only the claimant counts as auxiliary series. In general, the models with Google Trends tend to achieve a better estimation and nowcast of the change in unemployment, $R^{k,y}_t$, rather than the other two state variables, with respect to the models that include the claimant counts.
Including two instead of one factor clearly increases the complexity of the model, which is reflected in smaller accuracy gains (in the CC & GT model probably also due to the decreased magnitude of the correlation parameter with the claimant counts), especially for the nowcast of the state variables, with respect to including only one factor. Nonetheless, the correlations with both factors are individually and jointly significantly different from zero, indicating that both factors bring additional information about the Dutch unemployment.
Notice that in general all the relative measures of accuracy are below one, indicating that both the claimant counts and the Google Trends improve the estimation and nowcast accuracy of the unemployment and its change. Even when the Google Trends are not targeted and their factor is not significantly related to the unemployment, the measures are never drastically above one, meaning that our method tends to ignore auxiliary series that are not related to the target variable.
Finally, when we specified the covariance matrix (ref) in Section (ref), we did not let the claimant counts and the Google Trends be correlated because our goal is to improve the estimation/nowcast accuracy of the unobserved components of the labour force series, not of the claimant counts nor the Google Trends. Nonetheless, if the state variables of equation (ref) are all cointegrated (i.e. the correlation parameters are all equal to one) a more efficient estimation method would be to only estimate the variance of their common source of error. We therefore estimate the CC & GT model with one factor, when all series are correlated. We call this model “CC & GT all corr.". Table (ref) reports the empirical results also for this model. Although the nowcast accuracy is similar to the same model without the additional correlation between the claimant counts and the Google Trends (which we indicate as $\rho_{1,CC,GT}$), the in-sample accuracy deteriorates (even with respect to the baseline model), and $\rho_{1,CC,GT}$ is not significantly different from zero. We therefore conclude that the specification of the covariance matrix (ref) is appropriate.
In Table (ref) we report the empirical results for the GT and CC & GT models which employ the targeted Google Trends observed at the weekly frequency, and aggregated to the monthly frequency according to equation (ref) in order to include them in the models. In this case we still look at the sensitivity of the results with respect to the number of factors included in the model, but also with respect to the two additional methods for the estimation of $\bm \varLambda$ and $\bm \varPsi$ discussed at the beginning of this section.
The measures of accuracy are again broadly lower than one, but the gains are not as large as observed for the monthly Google Trends. Including two factors improves the accuracy in the GT model, but not in the CC & GT model, except for a more precise nowcast of $R^{k,y}_t$. The correlation parameter with the claimant counts remains large and significant. On the contrary, the correlation parameter with the first factor of the Google Trends is not significantly different from zero, and there is a weak evidence for the second factor being significantly related to the change in unemployment. For this reason we continue the analysis by considering two factors in the model.
Estimating $\bm \varLambda$ and $\bm \varPsi$ on the weekly Google Trends improves the measures of accuracy only for the CC & GT model, and not for the GT model. An additional iteration of the two step estimator, in order to obtain more accurate estimates of $\bm \varLambda$ and $\bm \varPsi$, achieves instead better nowcasts for both the GT and the CC & GT models (and also better in-sample estimates for the latter model), and a similar performance to the models which employ the monthly Google trends and include two factors. Notice that the values of the log-likelihood for these two models increased with respect to the same model specifications that use the original two-step estimation (without the additional iteration). The latter result, as pointed out in the explanation of the iterated estimation of $\bm \varLambda$ and $\bm \varPsi$ at the beginning of this section, is to be expected. Despite the above-mentioned improvements in estimation/nowcast accuracy, the correlation parameters with the Google Trends' factors are always insignificant. The aggregation of the Google Trends from the weekly to the monthly frequency yields time series that are more noisy with respect to the Google Trends that are directly observed at the monthly frequency, and detecting significant results therefore becomes harder.
Finally, even though weekly Google Trends allow to perform the monthly nowcasts on a weekly basis, we notice that, in general, the precision of the nowcast does not monotonically improve with the number of weeks. If the high-dimensional state space model could be expressed and estimated on the highest frequency, the weekly gains in nowcast accuracy could be more evident. Nonetheless, we are limited by the transition equations for the RGB and the survey errors, to estimate the model on the monthly frequency.
Figures (ref)-(ref) compare the point nowcasts, respectively, of the change in unemployment, its trend, and the population parameter, obtained with the baseline, the CC, and the GT and CC & GT models which employ monthly Google Trends and include two of their factors. From the first graph, it is evident that the models including claimant counts tend to deviate from the baseline model. The latter, on the contrary, gives similar results as those of the GT model. The point nowcasts of $L^{k,y}_t$ and $\theta^{k,y}_t$ are more similar throughout the model specifications, with a slight and positive difference between the models that include the Google Trends and the ones that do not, at the beginning of the out-of-sample period.
Figures (ref) and (ref) show the selection frequency of, respectively, the monthly and weekly Google Trends in the out-of-sample period. Some of the most selected search terms in both cases are: werklozen (unemployed people), baan zoeken (job search), curriculum vitae voorbeeld (curriculum vitae example), ww uitkering (unemployment benefits), ww aanvragen (to request unemployment benefits), resume, tijdelijk werk (temporary job), huizenmarkt zeepbel (housing market bubble). Notice that the latter term (as well as “economische crisis" (economic crisis) or “failliet" (bankrupt), which are also frequently selected monthly Google Trends) is of economic uncertainty nature, rather than being job-search related. A previous version of this paper only used the latter type of search terms, and did not find them to have explanatory power for the Dutch unemployment, which is now instead significantly improved by the inclusion of search terms related to economic uncertainty.
The results of the empirical analysis can be summarized as follows. Targeting the Google Trends improves the explanatory power of the latter series for the Dutch unemployment. Monthly Google Trends significantly improve the estimation and nowcast accuracy of the Dutch unemployment and its change, with both one and two factors. The largest gains are obtained when both the claimant counts and the Google Trends are included, and considering only one factor for the latter series. When two factors are considered, the gains are smaller but both factors seem to be significantly related to the change in unemployment, indicating that both of them should be included in the model in order to exploit all the information that the Google Trends give about the target variable. The sensitivity to the number of factors is somewhat similar for the weekly Google Trends, although there is a weak evidence only for their second factor to have a significant relation with the change in unemployment. The weekly Google Trends are less informative about the Dutch unemployment, yielding in general less improvements in estimation and nowcast accuracy, with respect to the monthly Google Trends. The contributions of the two types of Google Trends are comparable only when the two-step estimator is additionally re-iterated for the weekly Google Trends (in order to obtain more precise estimates of $\bm \varLambda$ and $\bm \varPsi$). This result suggests that iterating the two-step estimation can improve the explanatory power of the Google Trends, and that the latter series are sensitive to the estimates of $\bm \varLambda$ and $\bm \varPsi$. Improvements are, instead, not always present when $\bm \varLambda$ and $\bm \varPsi$ are estimated on the weekly data. In general, the claimant counts mainly have a positive impact on the estimation and nowcast accuracy of $\theta^{k,y}_t$ and $L^{k,y}_t$, whereas the Google Trends on $R^{k,y}_t$. The point nowcasts of the latter state variable are more sensitive to the type of auxiliary series included, with respect to the ones of $\theta^{k,y}_t$ and $L^{k,y}_t$.
The assumptions of normality made and discussed throughout the paper can be tested on the standardized one-step ahead forecast errors durbinkoopman2012: $\tilde{\bm v}^k_t = \bm B^k_t \bm v^k_t$, for $t=d+1, \dots, T$ with $(\bm F^k_t)^{-1}=\bm B^{k'}_t \bm B^k_t$, where $\bm F^k_t$ is the covariance matrix of the prediction errors $\bm v^k_t$ estimated with the Kalman filter. The prediction errors for the labour force are defined as $\bm v^{k,y}_t = \bm y^k_t - \bm Z^y_t \hat{\bm \alpha}^{k,y}_{t|\Omega_{t-1}}$, for the claimant counts as $v^{k,CC}_t = x^{k,CC}_t - \bm Z^{CC} \hat{\bm \alpha}^{k,CC}_{t|\Omega_{t-1}}$, and for the Google Trends as $\bm v^{k,y}_t = \bm x^{k,GT}_t - \hat{\bm \varLambda} \hat{\bm f}^{k}_{t|\Omega_{t-1}}$, for $t=d+1, \dots, T$ (the expressions for $\bm Z^y_t$ and $\bm Z^{CC}$ can be found in Appendix (ref)). We test the assumptions on the estimated CC & GT models when two factors of the Google Trends are included, and which employ, respectively, the monthly Google Trends, and the weekly Google Trends with the additional iteration of the two-step estimator (as they yield the best results in terms of estimation and nowcast accuracy of the state variables of interest, when two factors of the Google Trends are included).
We test the null hypothesis of univariate normality for each of the prediction error, with the ShapiroWilk1965 and BowmanShenton1975 tests, as suggested, respectively, in harvey_1990 and durbinkoopman2012. The former test is based on the correlation between given observations and associated normal scores, whereas the latter test is based on the measures of skewness and kurtosis.
The p-values from the Shapiro-Wilk test are reported in Figures (ref) and (ref) for the two different model specifications discussed above, respectively. For both model specifications, there is no (strong) evidence against the normality assumptions for the error terms of the labour force and the claimant counts series, as their corresponding p-values are above the confidence level of 0.05. This result suggests that the model is correctly specified for these series. The test instead rejects the null hypothesis of normality for most of the idiosyncratic components of the Google Trends. The normality assumption seems therefore not appropriate for the latter series, but as discussed in Sections (ref) and (ref), and examined in the simulation study of Appendix (ref), this type of misspecification does not affect the consistency of the estimators of the state variables and the hyperparameters, and does not seem to influence the performance of our method, nor the distribution of the LR test which allows to perform inference on the correlation parameters\footnote{Notice that we do not control for multiple hypotheses testing in this case. If we would control for it, we would obtain less rejections of the null hypothesis of normality for the error terms of the Google Trends, but the conclusions for the error terms of the labour force and the claimant counts series would stay the same.}. The conclusions from the Bowman-Shenton test are the same and the corresponding p-values are reported in Figures (ref) and (ref).
This paper proposes a method to include a high-dimensional auxiliary series in a state space model in order to improve the estimation and nowcast of unobserved components. The method is based on a combination of PCA and Kalman filter estimation to reduce the dimensionality of the auxiliary series, originally proposed by Doz2011, while the auxiliary information is included in the state space model as in harvey2000. In this way we extend the state space model used by Statistics Netherlands to estimate the Dutch unemployment, which is based on monthly LFS data, by including the auxiliary series of claimant counts and Google Trends related to job-search and economic uncertainty. The strong explanatory power of the former series, in similar settings, has already been discovered in the literature (see harvey2000 and Brakel2016). We explore to which extent a similar success can be obtained from online job-search and economic uncertainty behaviour. The advantage of Google Trends is that they are freely available at higher frequencies than the labour force survey and the claimant counts, and, contrary to the latter, they are not affected by publications delays. This feature can play a key role in the nowcast of the unemployment, as being the only real-time available information.
A Monte Carlo simulation study shows that in a smooth trend model our proposed method can improve the $\operatorname{MSFE}$ of the nowcasts of the trend's level and slope up to, respectively, around 25% and 75%. These results are robust to misspecifications regarding the distribution of the idiosyncratic components of the auxiliary series. Therefore, our method does have the potential to improve the nowcasts of unobserved components of interest.
In the empirical application of our method to Dutch unemployment estimation and nowcasting, we find that our considered Google Trends (when first targeted with the elastic net) do in general yield gains in the estimation and nowcast accuracy (respectively up to 40% and 25%) of the state variables of interest, with respect to the model which does not include any auxiliary series. This result stresses the advantage of using the high-dimensional auxiliary series of Google Trends, despite involving a more complex model to estimate, which is especially relevant for countries that do not have any data sources related to the unemployment (such as the registry-sourced series of claimant counts), other than the labour force survey. We also find that, under certain model specifications, including both claimant counts and Google Trends outperforms the model which only includes the former auxiliary series. This result is explained by the fact that the two auxiliary series have a positive impact on the estimation/nowcast accuracy of different unobserved components which constitute the unemployment, thus yielding an overall improvement of the fit of the model. This also indicates that claimant counts and Google Trends do not bring redundant information about the Dutch unemployment.
The magnitude of the above-mentioned gains is, nonetheless, sensitive with respect to the following aspects of the data and the model specification. First, in our empirical application we employ both monthly and weekly Google Trends. The latter need to be aggregated to the monthly frequency in order to be included in the model, but allow to perform the nowcast on a weekly basis. We find that the former are less noisy and provide in general more accurate estimates/nowcasts of the state variables of interest, with respect to the latter. The explanatory power of the monthly Google Trends for the Dutch unemployment is further corroborated by results from LR testing, which are in favour of their inclusion in the model. There is, instead, not strong and consistent evidence for this when the weekly Google Trends are employed.
Second, PCA involves the estimation of common factors that drive the Google Trends, and in our method we relate these factors to the unobserved components that constitute the Dutch unemployment. Information criteria suggest that the Google Trends are driven by either one or two common factors. We find that including two factors yields, in general, less gains in accuracy, with respect to including one factor (due to the increased complexity of the model), but there is evidence that the second factor is related to the unemployment, and therefore it should be included in the model in order to exploit all the information that the Google Trends give about the unemployment.
Finally, our estimation method is based on a two-step procedure. In the first step, the matrix of factors' loadings and the covariance matrix of the idiosyncratic components of the Google Trends are estimated by PCA. In the second step, these matrices are replaced by their PCA estimates, in order to re-estimate the Google Trends' factors and the unobserved components of the labour force series, with the Kalman filter. Replacing these matrices by their estimates might affect the explanatory power of the Google Trends. We find that the explanatory power of the weekly Google Trends can be improved (in order to yield similar gains as the ones obtained with the monthly Google Trends), with an additional iteration of the two-step estimation procedure, which should provide more accurate estimates of the two matrices.
As already mentioned, we generally find estimation/nowcast accuracy gains from the inclusion of the Google Trends, when they are first “targeted", by selecting the ones that are relevant for the Dutch unemployment, based on the elastic net penalized regression. If the targeting is not first applied, we do not find gains and significant relationships between the Google Trends and the Dutch unemployment. Nonetheless, in this case the results do not deteriorate with respect to the model that does not include any auxiliary series, suggesting that our method is able to ignore the inclusion of irrelevant auxiliary series, in the estimation/nowcast of unobserved components of interest. This result is corroborated in our Monte Carlo simulation study. Hence, our proposed approach provides a framework to analyse the usefulness of “Big Data" sources, with little risk in case the series do not appear to be useful.
One limitation of the current paper is that it does not allow for time-variation in the relation between the unobserved component of interest and the auxiliary series. For example, legislative changes may change the correlation between unemployment and administrative series such as claimant counts. Additionally, one can easily imagine the relevance of both specific search terms as well as internet search behaviour overall to change over time. While such time-variation may partly be addressed by considering shorter time periods, decreasing the already limited time dimension will have a strong detrimental effect on the quality of the estimators. Therefore, a more structural method is required that extends the current approach by building the potential for time variation into the estimation method directly, while retaining the possibility to use the full sample size. Such extensions are currently under investigation by the authors.