1. Introduction

In recent years, national statistical institutes have increasingly adopted non-traditional data sources for the production of official statistics. These data sources are not originally intended for statistical purposes, as they are collected for objectives other than the generation of official statistics. For instance, point-of-sale (POS) data recorded in retail stores, real-time population data derived from mobile phone GPS logs, price information aggregated from hundreds of online retailers, satellite imagery, and internet search trends—such as those published in Google Trends—can all be utilized to construct official statistics.

These emerging data sources have the potential to significantly enhance the utility of official statistics. If economic conditions can be captured more comprehensively and in a timelier manner through such sources, it is reasonable to expect improvements in the quality of decision-making, including economic policy formulation.

In the era of big data, private companies accumulate vast volumes of transaction data. It is therefore imperative to explore avenues for incorporating these data into the production of official statistics. However, the pace of progress in utilizing such data has been slower than anticipated. This is primarily because most of these data cannot be directly employed for statistical purposes. As they are not collected through survey methodologies grounded in sample design, they exhibit substantial selection bias when considered as potential sources for official statistics. In other words, if this pronounced selection bias can be appropriately corrected, a considerable volume of valuable transaction data could be effectively leveraged as a basis for official statistics, thereby expanding their scope and enhancing the quality of diverse policy and business decisions.

Previous studies ^1,2 have proposed frameworks for producing official statistics using new data sources subject to selection bias—particularly transaction data collected for operational rather than statistical purposes and not derived from survey-based methodologies. Specifically, prior work introduced a framework for generating preliminary reports based on partial datasets that, while biased, offer high timeliness. This approach is especially relevant for indicators where comprehensive nationwide data are eventually available for official statistics, but long publication lags reduce their usefulness for timely policy responses. In this context, the previous study emphasizes timeliness as a key advantage of new data sources. Surveys designed for the production of official statistics are inherently time-consuming due to the need to collect responses from firms or households. In contrast, the proposed framework can access new data sources with minimal lag—yesterday’s data become available the next day. If this framework enables the production of official statistics that reflect current economic conditions more promptly, it is reasonable to anticipate improvements in the quality of decision-making, including economic policy planning.

The problem of inferring population-level characteristics from biased data samples is a well-known challenge in machine learning, and extensive research has addressed this issue. In this light, integrating machine learning techniques into the production of official statistics holds the potential to fundamentally resolve longstanding challenges in the field. The aforementioned study focuses on methodologies rooted in machine learning—specifically, density ratio estimation and the concept of covariate shift.

In the present study, we demonstrate that the application of deep learning to the approach proposed in ^1,2 can further enhance estimation accuracy. In doing so, we adopt the methodology introduced in ³, which formulates a framework for Bregman divergence minimization. This framework offers a unified view of various density ratio estimation methods and investigates its applicability to highly flexible models such as deep neural networks. According to ⁴, recent advances in machine learning have highlighted the superior performance of deep neural networks in diverse domains including computer vision ⁵ and natural language processing ⁶, thus motivating their application to density ratio estimation. However, employing Bregman divergence minimization with highly flexible models introduces the risk of overfitting, particularly in the form of train loss hacking, which arises from properties of empirical Bregman divergence estimators. To mitigate this issue, ⁴ proposes a non-negative correction to the empirical estimates, thereby stabilizing the training process and reducing overfitting.

Following the methodology of ^1,2, we experimentally confirm that prompt preliminary statistics can be generated using biased data from a Japanese private employment agency.

For policymakers, corporate decision-makers, and recruitment professionals, timely insights into labor market trends are essential. To support this need, various official labor statistics are available in many countries, offering metrics such as unemployment rates, labor force counts, and average wages—typically published with a lag of two to three months. However, some indicators remain difficult to capture in a timely manner, with publication delays exceeding one year. One such indicator is the wage changes upon job change, capturing the pressure exerted on wages by supply–demand conditions in the external labor market. The timely release of this indicator would greatly benefit those involved in public policy, corporate management, and recruitment.

Among the various potential measures of the wage changes upon job change, this study adopts the proportion of job switchers who experience a wage increase, specifically those whose post-transition wages rise by more than 10%. The indicator is calculated as the ratio of individuals with a wage increase exceeding 10% to the total number of job switchers within a given period. The “more than 10%” threshold signifies a significant increase, and thus the indicator reflects the proportion of individuals whose wages markedly improved following a career change. In official Japanese labor statistics, this indicator is released with a time lag ranging from 6 to 13 months. The central challenge addressed in this study is the elimination of this extended delay.

Career transitions in the external labor market occur through multiple channels, including job advertisements, public and private employment agencies, and personal referrals. In this study, we focus on transaction data maintained by private employment agencies, as these organizations possess such data in formats suitable for real-time analysis. Private agencies typically conduct interviews with job seekers, enabling accurate documentation of pre-transition wages. Moreover, as part of the job change process, these agencies calculate and present to job seekers the annual wages they would earn in the new position. Consequently, they inherently possess both pre- and post-transition wage data. In contrast, job transitions occurring via advertisements or referrals generally lack mechanisms for recording wage data. While public employment agencies, in principle, could collect such information for statistical purposes, this system has not been fully implemented in Japan.

We obtained transaction data from the largest private employment agency in Japan. Additionally, we were granted access to anonymized sample data from official government statistics—unavailable to the general public—for benchmarking purposes. Both datasets consist of raw micro-level data rather than pre-tabulated statistics. Most of the processing in this study is conducted at the individual data level, not on aggregated statistical tables. It is important to emphasize that individuals cannot be identified and that strict data protection protocols have been observed throughout this study.

Although transaction data from private employment agencies are available in real time, their population coverage is limited—less than 5% even for the largest agency used in this study. Moreover, when used as a sample to infer trends at the national level, these data exhibit considerable bias. For example, in the second half of 2018, the average age of full-time job changers in the official statistics was approximately 40 years, whereas it was around 31 years in the transaction data. Similarly, the share of individuals with a university degree or higher was about 35% in the official data, compared to 79% in the transaction data. These discrepancies highlight the substantial sampling bias inherent in the transaction data.

This bias arises because the sample is limited to individuals who changed jobs through private employment agencies—namely, those firms who are willing to hire even when paying agency fees. In contrast, although official statistics are subject to significant publication delays, they capture a comprehensive view of national labor market dynamics. Therefore, if the relationship between private agency transaction data and official statistics can be learned and used to correct the bias in the former, the challenge can be effectively addressed.

The previous studies ^1,2 demonstrated the feasibility of correcting selection bias through density ratio estimation combined with covariate shift correction. Building upon that foundation, the present study introduces deep learning to further improve estimation accuracy. The main contributions of this study are as follows.

Although deep learning is generally considered prone to overfitting when applied to density ratio estimation, this study improves estimation accuracy without inducing overfitting.

It further shows that deep learning can be applied effectively without extensive hyperparameter tuning, demonstrating its practical applicability to the production of official statistics.

It should be noted that the proposed framework is not limited to labor statistics, but can also be applied to other statistical domains that exhibit similar structural biases.

2. Related Work

As noted in the previous section, national statistical institutes have recently begun adopting non-traditional data sources for the production of official statistics. These new data sources are not inherently designed for statistical purposes. The transaction data from private employment agencies used in this study fall into this category. Other examples have also been documented in the literature, where a variety of alternative data—such as satellite imagery, GPS-based population data from mobile phones, POS records from retail outlets, online price information, and internet search trends like Google Trends—have been employed in the production of official statistics.

The International Telecommunication Union (ITU), a United Nations agency specializing in information and communication technology, has been at the forefront of utilizing GPS-based real-time population data from mobile phones. As the secretariat of the Mobile Data Task Team within the United Nations Statistics Division, the ITU has promoted the global adoption of such data sources ⁷. Between 2016 and 2018, countries including Colombia, Georgia, Kenya, the Philippines, Sweden, and the UAE collected and processed data from telecom providers to complement household surveys and official registries, yielding 16 distinct indicators. During 2020–2021, Brazil and Indonesia conducted studies leveraging telecom data to assess two Sustainable Development Goal indicators: mobile network coverage and internet access. To support broader adoption, the ITU published the Handbook on the Use of Mobile Phone Data for Official Statistics ⁸. In Estonia, due to budget constraints, traditional travel surveys were replaced with mobile SIM card data to monitor travel expenditures ⁹.

Since the early 2000s, the use of POS data has expanded significantly across numerous countries. Notably, Switzerland, Norway, and the Netherlands have undertaken dedicated efforts to incorporate POS data as a way to overcome limitations of conventional data sources. Beginning with the study cited in ¹⁰, new price indices based on POS data have been developed and validated in subsequent research ^11,12. In Japan, a private company provides price indices constructed from POS data in accordance with the methodology described in ¹³. These indices are utilized by the Bank of Japan in various official reports. Additionally, price indices have been created using price data collected from a wide range of online retailers. The Billion Prices Project, launched in 2008 by Professors Alberto Cavallo and Roberto Rigobon at MIT Sloan and Harvard Business School, introduced novel price indices based on online retail data ^14,15, as well as a new methodology for constructing purchasing power parities (PPPs) ¹⁶.

The Organisation for Economic Co-operation and Development (OECD) has employed internet search data from Google Trends in its statistical initiatives. These data contribute to the publication of the OECD Weekly Tracker of Economic Activity, which estimates weekly GDP for 46 countries, including OECD and G20 members ¹⁷. Google Trends data have also been used in estimating unemployment rates ¹⁸, serving as an auxiliary time series for the nowcasting of monthly unemployment statistics alongside the official labor force survey.

The studies most closely related to the present work are ¹ and ², also discussed in Section 1. Among them, ² shares not only the methodological framework but also the experimental setting with the present study, and thus is used as a benchmark reference for comparison with this research. Both studies, however, are based on the same overall methodological framework. Those studies propose a framework for producing timely indicators of wage changes among job changers in Japan by utilizing transaction data from a private recruitment agency. The data sample is limited to individuals who changed jobs via private recruitment services, excluding those who transitioned through job advertisements, personal referrals, or public employment services. Compared to the datasets employed in the other cases described above, this results in a stronger degree of selection bias. To address this issue, ^1,2 apply methodologies developed in the field of machine learning, particularly density ratio estimation and the concept of covariate shift.

3. Task Setting

In this section, we describe the task structure, evaluation method, and the positioning of the present experiment within the typology of selection bias, following ². To maintain generality in the discussion, we do not elaborate on specific experimental settings or data details here; those will be provided in Section 5, following the methodological exposition in Section 4.

3.1. Structure of the Task

Official statistics are typically published in multiple stages, such as preliminary, revised, and final releases. Naturally, earlier releases tend to be less accurate, whereas later ones offer greater precision. In most cases, only the preliminary indicators are available in time to support economic policy and business decision-making. This structure is illustrated in Fig. 1, where the term error is used schematically to represent the discrepancy between the true statistical value—which is fundamentally unobservable—and the value obtained through surveys.

However, many important indicators in official statistics do not include preliminary releases suitable for timely decision-making. This is often due to difficulties such as the high cost of surveys or the burden placed on businesses, which complicate rapid data collection and publication. The present study focuses on such indicators. As shown in Fig. 1, for the indicator used in our experiment, only the final report is publicly available. In this study, we demonstrate that it is possible to construct a preliminary report even for such indicators.

In Fig. 1, the term error schematically represents the discrepancy between the true, unobservable statistical value and the survey-based observed value. However, because the true values are unobservable in this study, we treat the final released statistics as ground truth and evaluate the accuracy of our preliminary estimates by measuring the difference between the estimated and final values. To assess estimation performance, we employ the mean absolute error (MAE). Specifically, the MAE is calculated by comparing our estimated proportion of individuals who experienced wage increases after changing careers with the corresponding values published by the government.

As mentioned earlier, this study aims to improve estimation accuracy by incorporating deep learning into the framework proposed in ². Accordingly, we evaluate whether the MAE is further reduced under the same experimental settings used in the previous study.

3.2. Evaluation

We focus on MAE because it provides an interpretable measure of average deviation in percentage points and is robust to extreme observations, which is particularly relevant in macroeconomic time-series settings with a limited number of evaluation periods. For reference, we also confirmed that the two key findings reported in Section 6 are qualitatively similar when alternative evaluation metrics such as MSE and RMSE are used (results omitted for brevity).

Both the previous and present studies consider the following estimation strategies: (a) using the estimated density ratio solely as a weight, (b) incorporating an additional supervised learning step (classification), and (c) incorporating an additional supervised learning step (regression). In Fig. 3 in Section 4, the leftmost panel corresponds to (a), the center panel corresponds to (b), and the right panel corresponds to (c). The comparative framework underlying these strategies is illustrated in Fig. 2.

As in the prior work, we also report a simple extrapolation-based prediction as a baseline, which is defined by the following equation. This baseline reflects a setting in which no correction for selection bias is applied:

\begin{equation*} \hat{o}_{T} = o_{T-2} \dfrac{s_{T}}{s_{T-2}}, \end{equation*}

where \(\hat{o}_{T}\) denotes the predicted value of the official statistics for period \(T\), and \(o_{T-2}\) is the actual value from period \(T-2\). \(s_{T}\) and \(s_{T-2}\) are supplementary indicators for periods \(T\) and \(T-2\), respectively. In this study’s empirical setting, both \(o\) and \(s\) represent the proportion of individuals whose wages increased by 10% or more after changing jobs.

In relation to the notation used in Section 4, \(o_{T}\) corresponds to the number of components equal to 1 in the vector \(Y_{g,T}\), divided by the number of dimensions. Likewise, \(s_{T}\) corresponds to the number of components equal to 1 in the vector \(Y_{p,T}\) (under a classification label setting), divided by the number of dimensions.

4. Methodology

This section presents the methodological framework of our study. Fig. 3 illustrates the overall architecture of the proposed estimation procedure. This study proposes the application of deep learning to the step corresponding to the second layer from the top in Fig. 3. Section 4.2 outlines the methodological framework common to the prior studies ^1,2. Subsequently, Section 4.3 provides a detailed description of the newly introduced deep learning approach adopted in this study. Accordingly, Section 4.3 represents the primary methodological contribution and the key difference from the existing studies. Note that the prior studies ^1,2 provide more detailed explanations with figures for Sections 4.1, 4.2, 4.4, 4.5, and 4.6; therefore, we only provide a brief overview here.

4.3. Incorporating Deep Neural Networks into Density Ratio Estimation

In the previous studies ^1,2, uLSIF ²⁵ was adopted for density ratio estimation, where the model employed was linear-in-parameter, with nonlinear structures captured through kernel functions and regularization introduced to prevent overfitting. In contrast, the model in this study is a deep neural network, and we investigate the potential for improving the accuracy of density ratio estimation. Section 4.3.1 presents the concept of density ratio matching based on Bregman divergence, showing that uLSIF is a special case within this framework. Section 4.3.3 discusses the overfitting issue known as train loss hacking, which arises when performing density ratio matching with deep neural networks under the Bregman divergence framework. Finally, Section 4.3.4 introduces the non-negative Bregman divergence estimator, a recently proposed method for improving the robustness of empirical Bregman divergence estimators against train loss hacking.

4.3.1. Density Ratio Matching by Bregman Divergence Minimization

Following ³, this section introduces the concept of density ratio matching using Bregman divergence, and shows that the uLSIF method presented in the previous section is a special case within this broader framework. This formulation generalizes the least-squares approach and also encompasses other methods examined in ^1,2, such as moment matching ^26,27, probabilistic classification ^28,29,30, density matching ^{31,32,33,34,35}, and density ratio fitting ²⁵.

The objective is to estimate the density ratio:

\begin{equation} \label{eq:dens} w(x) = \frac{p_T(x)}{p_S(x)}, \end{equation}

from samples drawn from \(p_S(x)\) and \(p_T(x)\).

The core idea of density ratio matching is to directly fit a density ratio model \(\hat{w}(x)\) to the true density ratio function \(w(x)\) under a specified divergence. While this may appear similar to regression, density ratio matching differs fundamentally in that samples of the true density ratio are not available. In this work, we use the Bregman divergence to quantify the discrepancy between \(w(x)\) and \(\hat{w}(x)\).

The Bregman divergence is a generalization of Euclidean distance. Let \(f\) be a differentiable, strictly convex function. Then the Bregman divergence from \(u^*\) to \(u\) is defined as:

\begin{equation} \text{BR}'_f\left(u^*\|u\right) := f\left(u^*\right) - f(u) - \partial f(u)\left(u^* - u\right), \end{equation}

where \(\partial f\) denotes the derivative of \(f\).

The term \(f(u) + \partial f(u)(u^* - u)\) corresponds to the first-order Taylor expansion of \(f\) at \(u\), evaluated at \(u^*\). Thus, the Bregman divergence measures the difference between the true value of \(f(u^*)\) and its first-order linear approximation from \(u\). See Fig. 4 for illustration.

The discrepancy between the true density ratio function \(w\) and the model \(\hat{w}\) is measured using the Bregman divergence, defined as

\begin{align} \label{eq:63} \mathrm{BR}'_f \big(w \| \hat{w}\big) := \int & p_{S}(x) \biggl( f\big(w(x)\big) - f\big(\hat{w}(x)\big) \nonumber\\ & - \partial f\big(\hat{w}(x)\big)\big(w(x) - \hat{w}(x)\big)\biggr) dx. \end{align}

A key motivation for adopting this formulation is that the Bregman divergence enables a direct empirical approximation for any choice of \(f\). Specifically, we isolate the relevant component of \(\mathrm{BR}_f(\hat{w})\) as follows:

\begin{equation} \mathrm{BR}'_f \big(w \| \hat{w}\big) = \mathrm{BR}_f\big(\hat{w} \big) + C, \end{equation}

where

\begin{equation*} C := \int p_{S}(x) f\big(w(x)\big) dx, \end{equation*}

is a constant independent of \(\hat{w}\), and

\begin{align} \label{eq:65} \mathrm{BR}_f\big(\hat{w}\big) := \int p_{S}(x) \biggl( &\, \partial f\big(\hat{w}(x)\big)\, \hat{w}(x) - f\big(\hat{w}(x)\big) \biggr) dx \nonumber\\ & - \int p_{T}(x) \partial f\big(\hat{w}(x)\big) dx. \end{align}

An empirical approximation \(\widehat{\mathrm{BR}}_f(\hat{w})\) of \(\mathrm{BR}_f(\hat{w})\) is then given by

\begin{align} \label{eq:66} \widehat{\mathrm{BR}}_f\big(\hat{w}\big) :=\frac{1}{n} \sum_{i=1}^{n} \biggl(& \partial f\left(\hat{w}\big(x_{i}\big)\right) \hat{w}\left(x_{i}\right) - f\left(\hat{w}\left(x_{i}\right)\right) \biggr) \nonumber\\ &- \frac{1}{n'} \sum_{j=1}^{n'} \partial f\left(\hat{w}\big(x'_{j}\big)\right). \end{align}

This leads directly to the following optimization criterion:

\begin{equation} \min_{\hat{w}} \widehat{\text{BR}}_f\big(\hat{w}\big), \end{equation}

where \(\hat{w}\) is optimized over a predefined class of functions.

4.3.2. uLSIF Within the Framework of Density Ratio Matching

We now demonstrate that LSIF/uLSIF, as used in ^1,2, constitutes a special case of density ratio matching ³. Specifically, there exists a Bregman divergence for which the density ratio matching objective reduces to that of LSIF/uLSIF.

When

\begin{equation} f(u) = \frac{1}{2} (u - 1)^2. \end{equation}

Equation 3 is reduced to the squared (SQ) distance:

\begin{equation} \text{SQ}\left(u^* \| u\right) := \frac{1}{2} \left(u^* - u\right)^2. \end{equation}

Following Eqs. 5 and 6, we denote the squared loss without the irrelevant constant term by \(\text{SQ}(\hat{w})\), and its empirical approximation by \(\widehat{\text{SQ}}(\hat{w})\), respectively:

\begin{align} \mathrm{SQ} \big(\hat{w}\big) : &= \frac{1}{2} \int p_{S}(x) \hat{w}(x)^2 dx - \int p_{T}(x) \hat{w}(x) dx,\\ \end{align}

\begin{align} \widehat{\text{SQ}} \big(\hat{w}\big) : &= \frac{1}{2 n} \displaystyle\sum_{i=1}^{n} \hat{w}\left(x_i\right)^2 - \frac{1}{n'} \displaystyle\sum_{j=1}^{n'} \hat{w}\left(x'_j\right). \end{align}

4.3.3. Train Loss Hacking Problem

We now discuss the overfitting problem known as train loss hacking, which arises when performing density ratio matching with deep neural networks under the Bregman divergence framework. Given the success of deep neural networks in fields such as computer vision ⁵ and natural language processing ⁶, their application to density ratio estimation represents a natural extension. However, naive application within the Bregman divergence framework can lead to substantial overfitting. For instance, ⁴ demonstrates that LSIF implemented with neural networks exhibits severe overfitting. One potential cause of this issue has been hypothesized to be train loss hacking—a phenomenon also discussed in a related context by ³⁶ in the setting of positive-unlabeled (PU) learning.

The density ratio estimation problem has already been formulated in Section 4.3.1. The objective in this setting is to estimate the density ratio defined in Eq. 1, based on samples drawn from the source and target distributions \(p_{S}(x)\) and \(p_{T}(x)\).

Train loss hacking refers to a phenomenon in which \(\hat{w}(x)\) takes on excessively large values for \(\{ x'_j \}_{j=1}^{n'}~\mathrm{i.i.d.}~\sim p_{T}(x)\).

The structure of the train loss hacking problem is illustrated in Fig. 5.

This issue arises when a model is trained to minimize an objective function that involves multiple distinct sample sets. Specifically, in the case of density ratio matching under the Bregman divergence framework, the objective function defined in Eq. 6 consists of two empirical averages. Among these, the second term

\begin{equation*} - \frac{1}{n'} \displaystyle\sum_{j=1}^{n'} \partial f\left(\hat{w}\big(x'_{j}\big)\right), \end{equation*}

can be minimized by encouraging the model to assign excessively large values to specific input points \(\{x'_j \}_{j=1}^{n'} ~ \mathrm{i.i.d.}~\sim p_{T}(x)\).

Since \(f\) is a convex function, its derivative \(\partial f\) is monotonically increasing; thus, the second term decreases as \(\hat{w}(x'_j)\) increases. Consequently, when no lower bound exists for this term, it frequently diverges to negative infinity in numerical computations. Even if a bounded model is used for \(\hat{w}\), numerical instability often leads to divergence. Moreover, even when a lower bound is imposed, \(\hat{w}(x'_j)\) still tends to assume the maximum possible value within its output range at the points \(\{ x'_j \}_{j=1}^{n'}\). Simply capping the model’s output, for example, by applying \(\max\{\cdot, B\}\) with some constant \(B > 0\), does not resolve the issue, as the model still converges to its largest possible value.

This is a significant issue because merely increasing the output at \(\{ x'_j\}_{j=1}^{n'}\) does not constitute a valid training criterion for density ratio estimation, and it results in an unreasonable density ratio estimator. The problem becomes particularly severe in highly flexible models. If the hypothesis class is severely restricted, a trade-off may be introduced by the remaining term of Eq. 6:

\begin{equation*} \frac{1}{n} \sum_{i=1}^{n} \biggl( \partial f\left(\hat{w}\big(x_{i}\big)\right) \hat{w}\big(x_{i}\big) - f\left(\hat{w}\big(x_{i}\big)\right) \biggr). \end{equation*}

However, in the case of highly flexible models, such as deep neural networks, the model can easily fit \(\{ x'_j \}_{j=1}^{n'}\) and \(\{ x_i \}_{i=1}^{n}\) separately.

4.3.4. Deep Learning for Density Ratio Estimation Using Non-Negative Risk Estimator

Naive application within the Bregman divergence framework leads to a significant overfitting issue, known as train loss hacking, as discussed in Section 4.3.3. This section introduces the non-negative Bregman divergence estimator ⁴, a method proposed to enhance the robustness of the empirical Bregman divergence estimator against train loss hacking.

To address this problem, ⁴ proposed a non-negative Bregman divergence estimator that modifies the empirical objective defined in Eq. 6 to improve robustness. According to ⁴, the method is inspired by ³⁶, which introduces a non-negative correction to the empirical risk in PU learning, based on the observation that a portion of the population risk is inherently non-negative. However, in the context of density ratio estimation, it is not straightforward to apply this idea, since it is unclear which part of the population risk, as defined in Eq. 5, is non-negative. To overcome this, an upper bound \(\overline{W}\) is assumed for the density ratio \(w\), allowing the identification of a non-negative component of the population risk, as expressed in Eq. 3. A non-negative correction is then applied to the empirical Bregman divergence estimator based on this decomposition. This strategy also generalizes the concept of non-negative PU learning ³⁶. To ensure that this method effectively mitigates train loss hacking, the following assumption must be satisfied:

Assumption 1

The density ratio \(w\) is bounded from above, i.e.,

\begin{equation*} \overline{W} = \sup_{x \in \mathcal{D}}w(x) < \infty. \end{equation*}

Then, we specify a constant \(C\) such that \(0 < C < {1}/{\overline{W}}\). Based on this constant, we impose the following assumption.

Assumption 2

\(\tilde{f}\) defined by

\begin{equation*} \partial f (t) = C \big(\partial f(t)t - f(t)\big)+\tilde f(t), \end{equation*}

is bounded from above.

Then, the Bregman divergence minimization objective given in Eq. 5 can be reformulated as

\begin{align} \label{eq:613} \mathrm{BR}_f \big(\hat{w}\big) := & \int p_{S}(x) \left( \ell_{1}\big(\hat{w}(x)\big) \right) dx \nonumber\\ & - C \int p_{T}(x) \left( \ell_{1}\big(\hat{w}(x)\big) \right) dx \nonumber\\ & + \int p_{T}(x) \left( \ell_{2}\big(\hat{w}(x)\big) \right) dx \nonumber\\ & - \big(1-C\big)A, \end{align}

where \(\ell_{1}\) and \(\ell_{2}\) are

\begin{equation*} \begin{aligned} \ell_{1}(t) &:= \partial f (t) t - f (t) + A, \quad \ell_{2}(t) := -\tilde f(t), \end{aligned} \end{equation*}

and \(A\) is a constant such that \(\ell_{1}(t) \geq 0\) for all \(t\).

The first two terms in Eq. 12 are non-negative, as both \(\ell_{1}\) and \(p_{S} - C p_{T}\) are non-negative under Assumption 1 and the condition \(0 < C < {1}/{\overline{W}}\).

Motivated by this observation, the following modified empirical risk is proposed.

\begin{align} \mathrm{nn}{\widehat{\mathrm{BR}}}_f(\hat{w}) := & \biggl( \frac{1}{n} \sum_{i=1}^{n} \ell_{1}\big(\hat{w}\big(x_{i}\big)\big) \nonumber\\ &\quad-C \frac{1}{n'} \sum_{j=1}^{n'} \ell_{1}\big(\hat{w}\big(x'_{j}\big)\big) \biggr)_+ \nonumber\\ & + \frac{1}{n'} \sum_{j=1}^{n'} \ell_{2}\biggl(\hat{w}\big(x'_{j}\big)\biggr), \end{align}

where \((\cdot)_+ := \max\{0, \cdot\}\).

4.4. Supervised Learning Under Covariate Shift: Classification

We now perform supervised classification using the weighted private employment agency samples from Section 4.2. We aim to estimate a function \(F\) such that \(Y_{p,T} = F(X_{p,T})\), using the weights \(w_T\). Using the estimated function \(F\), we compute \(F(X_{g,T})\) and use the resulting classifier score in the subsequent step.

Three classifiers were tested in prior studies: logistic regression with an elastic net penalty, random forest, and gradient-boosted decision trees. Due to minimal differences in accuracy, we adopt only logistic regression with elastic net regularization in this study. All attributes listed in Section 5.2 are used as explanatory variables in \(X_{p,t}\).

Standard learning algorithms, such as maximum likelihood estimation, become biased under covariate shift, as described in Section 3.3. This bias can be asymptotically corrected using density ratio weighting ²⁰.

5. Experiment

While Section 3 maintained a general perspective by providing an abstract overview of the task setting, this section presents the specific details of the experimental design that were not previously discussed. In particular, we describe the time lag structure of the official statistics used, the characteristics of the datasets, and the full details of the experimental procedure.

5.3. Details of Experimental Design

To evaluate the performance of our estimation, we employed the MAE, as described in Section 3.2. The dataset used spans from the first half of 2004 to the second half of 2018, and the validation period ranges from the second half of 2013 to the second half of 2018.

For each target half-year within the validation period, we utilized government survey samples up to one year prior to the target half-year and private employment agency samples from the target half-year.

For example, as outlined in Section 4.1, we used government survey samples from the first half of 2004 to the second half of 2013 to estimate the second half of 2014 using SARIMA. Subsequently, as discussed in Section 4.2, we used private employment agency samples from the second half of 2014 to perform density ratio estimation via uLSIF.

In Section 4.6, the label information was corrected using a constant value \(\beta\), which was calculated as the ratio of two quantities: the numerator being the actual proportion of individuals with increased wages after changing careers, and the denominator being the corresponding estimated value without this correction.

We considered two approaches for computing \(\beta\). \(\beta\) is calculated using the ratio between the actual proportion of samples where \(Y=1\) in each period and its estimated value when label information is not corrected. In the first approach, \(\beta\) was calculated as the average of the ratios, which are the estimation results from the first half of 2013 to the most recent half-year preceding the target period. This setting ensures that only information available at the time of estimation is used. In the second approach, \(\beta\) was calculated as the average of the ratios, which are all the estimation results from the first half of 2013 to the second half of 2018. This setting assumes access to future data that would not have been available at the time of estimation and is therefore infeasible in practice. Under the first approach, relatively more information is available for later periods in the test window, whereas information is more limited for earlier periods. In contrast, the second approach uses richer and more uniform information across all periods. Thus, the second approach can be interpreted as a reference case that provides an upper bound on achievable performance under the assumption of time-invariant the ratio. However, this interpretation is valid only if the ratio is approximately time-invariant. If the ratio is unstable over time, the results under the second approach could be worse than those under the first approach.

6. Results

As mentioned in Section 3.2, this study aims to demonstrate improved estimation accuracy by applying deep learning to the approach proposed in the prior work ². Accordingly, we show that the MAE is further reduced under the same experimental settings as used in the prior study.

Both the prior study and the present work consider two estimation strategies: one that uses the estimated density ratio solely as a weight, and another that incorporates an additional supervised learning step. In Fig. 3 in Section 4, the leftmost case corresponds to the former approach, while the center and right cases correspond to the latter.

Therefore, in comparing the prior work and the present study, we examine whether applying deep learning improves estimation accuracy in each of these settings: (a) using density ratio weights only, (b) incorporating supervised learning (classification), and (c) incorporating supervised learning (regression). This comparison is illustrated in Fig. 2 in Section 3.2.

Tables 2 and 3 present the results under different hyperparameter \(1/C\) settings. The definition and role of the hyperparameter is described in Section 4.3.4. The upper row shows the results of (a) using density ratio weights only, while the middle row presents the results of (b) incorporating supervised learning (classification), and the lower row presents the results of (c) incorporating supervised learning (regression).

The best scores from the prior study, as also noted in the table footnotes, are as follows: for case (a) using density ratio weights only, 1.82 in Table 2 and 1.62 in Table 3; for case (b) incorporating supervised learning (classification), 1.60 in Table 2 and 1.46 in Table 3; and for case (c) incorporating supervised learning (regression), 1.70 in Table 2 and 1.49 in Table 3. As a reference, the MAE of the simple extrapolation-based prediction is 2.55.

The first notable finding is that the method proposed in this study outperforms the approach tested in the prior work, in terms of MAE. In Table 2, where \(\beta\) is calculated from the results before the target period, increasing the hyperparameter \(1/C\) beyond a certain level yields performance that surpasses that of the previous study across all three settings: (a) using density ratio weights only, (b) incorporating supervised learning (classification), and (c) incorporating supervised learning (regression). Moreover, comparing settings (b) incorporating supervised learning (classification) and (c) incorporating supervised learning (regression), setting (b) achieves lower MAE and more consistently outperforms the previous study, indicating its superior performance. In Table 3, we obtain similar results for setting (b).

It is worth noting that incorporating supervised learning yields better accuracy than using density ratio weights only, and that calculating \(\beta\) from all results leads to better accuracy than calculating it from the results before the target period, consistent with the findings of the previous study.

The second key finding concerns the behavior of MAE in response to changes in hyperparameter values. When the hyperparameter \(1/C\) is set to an excessively small value, MAE increases significantly. As \(1/C\) gradually increases, MAE correspondingly decreases. However, beyond a certain threshold, MAE stabilizes. Focusing on Table 2, which represents a more realistic setting, results that outperform the previous study begin to appear once \(1/C\) exceeds a certain level, and the results become more stable when \(1/C \geq 16\). In particular, all results in Table 2 outperform the previous study when \(1/C \geq 16\). These findings suggest that, in the present application, setting \(1/C\) to approximately 16 or slightly larger is sufficient in practice. This indicates that meticulous searching of this hyperparameter at a high computational cost is unnecessary; instead, it is sufficient to set it to a reasonably large value with ample margin. This finding is of particular importance in the context of official statistical practice, where methods involving substantial computational costs are generally considered impractical to implement.

7. Discussion and Conclusion

The objective of this study, along with previous studies ^1,2, is to develop a framework for producing official statistics using new data sources that are subject to selection bias. In recent years, national statistical institutes have increasingly begun to incorporate non-traditional data sources—such as POS data and mobile phone GPS data—into the production of official statistics. These initiatives are expected to enhance the quality of decision-making, particularly in the context of economic policymaking. Within this broader trend, the use of transaction data accumulated by private companies for official statistical purposes represents a natural progression. However, the adoption of such data has been slower than anticipated.

This slow uptake is primarily attributable to the significant selection bias inherent in these data when used as inputs for official statistics. If this strong selection bias can be appropriately corrected, then large volumes of otherwise underutilized transaction data could be effectively leveraged for statistical production. Doing so would substantially broaden the scope of official statistics and enhance the quality of decision-making across various domains.

Using Japanese labor statistics as a case study, this study and previous work ^1,2 demonstrate that even data affected by substantial selection bias can be used to construct informative preliminary indicators. As shown in Section 5.1, the specific indicator examined in this study typically suffers from a publication lag exceeding one year. Although inherently valuable, this delay renders the indicator ineffective for real-time decision-making. However, by applying the methodology proposed herein and in related studies, it becomes feasible to generate a timely version of the indicator that can support more responsive, real-time policy and business decisions.

This study explored the integration of deep learning into the density ratio estimation process. While the accuracy gains from incorporating deep learning are not immediately apparent—due to the overfitting issue known as train loss hacking, as discussed in Section 4.3.3—our experimental results suggest the potential for a certain degree of improvement.

Compared to the methods proposed in the previous studies ^1,2, the deep learning-based approach introduces greater computational overhead during training. Consequently, there exists a trade-off between accuracy gains and computational cost, and the choice between the conventional methods and the deep learning-based approach must be carefully weighed in practical implementations. Nevertheless, our findings indicate that despite the increase in training cost, the deep learning-based method does not require extensive hyperparameter tuning. Rather than conducting computationally expensive searches, it suffices to assign a reasonably large value with an adequate margin. This implies that the additional computational burden associated with the proposed method may be tolerable within the context of official statistical practice.

The contribution of this study is not confined to the labor statistics used in the empirical analysis. The proposed approach is applicable not only to these statistics but also to a broad range of other datasets that share similar structural characteristics. Indeed, many practical scenarios are expected to exhibit analogous data-generating mechanisms.

To date, the application of machine learning and deep learning techniques in the production of official statistics has remained extremely limited. Nevertheless, methods developed in these fields offer considerable potential—particularly in the context of utilizing non-traditional data sources. If the findings of this study contribute to the broader adoption of machine learning methodologies in official statistics, it would represent a meaningful advancement in the field.