High-dimensional censored MIDAS logistic regression for corporate survival forecasting

High-dimensional censored MIDAS logistic regression for corporate survival forecasting

\onehalfspacing

{\it Keywords:} Corporate survival analysis; high-dimensional censored data; mixed-frequency data; logistic regression; sparse-group LASSO

Introduction

Regulators, lenders, and investors are increasingly focused on identifying vulnerable firms and developing accurate models to predict firm failures well in advance, as the ability to correctly predict such failures could result in a more resilient financial stability policy and better financial outcomes for market participants. As a result, an extensive body of literature is dedicated to understanding the determinants of firm failures. Traditional statistical models, such as discriminant analysis almon1965distributed, logistic regression ohlson1980financial, and hazards models shumway2001forecasting, along with other time-sensitive approaches duffie2007multi, have historically been the main focus of study. However, with the advent of more extensive datasets in recent years, the focus has increasingly shifted toward machine learning methods, which are better equipped to handle high-dimensional data. Such models have shown superior accuracy in predicting firm failures barboza2017machine due to their efficient handling of rich data sources. Over time, the task of forecasting corporate survival has gained significant attention due to its critical economic implications and its close connection to other challenges, such as predicting household loan defaults.

In this paper, we focus on the task of predicting the probability that a firm will fail within the first $t$ years after its initial listing, conditional on its survival for the first $s$ years, where $s < t$. This problem presents three significant statistical challenges. First, data are often right-censored. {\color{black} To illustrate how censoring arises in practice, we consider our empirical dataset on Chinese listed manufacturing firms (see Section (ref) for details). Figure (ref) depicts firms with different censoring statuses over the observation period from January$1^{\text{st}}$, 1985 to December $31^{\text{st}}$, 2020. For each firm, the survival time $T$ is defined as the number of years from its initial public offering (IPO) date to the first occurrence of financial distress. If a firm is never classified as distressed during the observation period, its survival time $T$ is right-censored. In this case, $T$ exceeds the censoring time $C$, and only $C$, defined as the number of years from the IPO date to the end of the observation period, is observed. In Figure (ref), Firms $1$ and $2$ are uncensored, so their survival times are fully observed, whereas Firms $3$ and $4$ are censored and only their censoring times $C$ are observed.

} Second, the high dimensionality of the predictors adds complexity. Modern data sets provide a wealth of variables for each listed firm, increasing the analytical burden. Third, the mixed-frequency nature of the data compounds the difficulty. For each potential predictor, we observe numerous lags, exacerbating the challenge of managing the proliferation of parameters.

As highlighted in our review of the literature below, in our view, the existing methods for this prediction task do not adequately address all three challenges. To bridge this gap, we propose a novel high-dimensional censored MIDAS logistic regression method that addresses these complexities. Our approach is based on a high-dimensional logistic regression framework to estimate survival probabilities. To address right-censoring, we make use of a tool from the survival analysis literature called outcome-weighted inverse probability of censoring weighting, as described in blanche2023logistic. The mixed-frequency nature of the data is managed using mixed data sampling (MIDAS), an approach developed and popularized by ghysels2007midas. This method approximates the coefficients of the lags of each variable using a finite-dimensional series basis, known as the dictionary. Finally, to handle the high dimensionality of the predictors, we apply a sparse-group LASSO penalty. This penalty not only manages the dimensionality of the regressors but also accounts for the group structure of the predictors, which corresponds to the lags of the original variables, as discussed in babii2022machine.

We derive finite-sample bounds on the estimation error of our penalized estimator. Notably, these bounds allow for heavy-tailed variables and account for both the approximation error and right-censoring, which are novel contributions to the literature on high-dimensional logistic regression models. \textcolor{black}{Additionally, we propose a de-sparsified procedure for the penalized estimator, which leverages the nodewise regression van_de_geer_debias. Under suitable regularity conditions, the de-sparsified estimator is asymptotically normal, enabling valid inference on individual components of the parameter vector while properly accounting for right-censoring, MIDAS approximation error, and heavy-tailed covariates. Notably, the presence of censoring introduces a nonstandard variance structure for the de-sparsified estimator, which has not yet been examined in the existing literature van_de_geer_debias,caner2023generalized,babii2023machine.}

The finite-sample performance of our method is evaluated through simulations, demonstrating its robustness against natural alternatives. Furthermore, we showcase the practical advantages of our approach through an application to forecasting the financial distress of Chinese listed firms. In this context, our method significantly outperforms the standard logistic regression benchmark and other competing methods over several horizons, underscoring its empirical effectiveness. Several practical augmented prediction methods, including oversampling and incorporating macro data into the model, are utilized. To further demonstrate the effectiveness of the proposed method that includes censoring information, a comparison is conducted with a method that excludes censored firms. \textcolor{black}{The de-sparsified estimator is also employed to identify statistically significant financial predictors of firm distress. These identified predictors provide practitioners with new empirical insights and help direct attention to specific covariates that are particularly informative for distress prediction.} Finally, our novel approach is implemented in the R package Survivalml to make it readily accessible for practitioners.\footnote{The package is publicly available at \url{https://github.com/Wei-M-Wei/Survivalml}.}

\paragraph{Literature review.} Let us first review how the existing methods address the three challenges we described, which are inherent in corporate survival forecasting. This review will stress the advantages of our methodology over popular alternatives. Given the extensive literature on this topic, we do not aim to provide an exhaustive review. Instead, we focus on surveying key approaches to the problems at hand. We also cite papers on the related problem of forecasting loan default.

To address the right-censoring of data, many studies restrict their analysis to firms that were first listed more than $t$ years before the end of the follow-up period audrino2019predicting,petropoulos2020predicting. Under the classical assumption of independent censoring, this approach avoids selection bias. However, it discards data on firms listed less than $t$ years ago, leading to a loss of efficiency. Another common strategy is to directly model the hazard rate of firm failure, using methods such as Cox models or single-index models ding2012class,lee2014business,kim2016survival,zhou2022recurrence,li2023corporate. Although effective in some contexts, this approach has limitations. Typically, the primary interest lies in estimating the probability of failure, not the hazard rate, and, therefore, modeling the survival probability as we do is more natural to solve the problem at hand. Furthermore, none of the aforementioned hazard-based approaches explicitly account for the challenges posed by high-dimensional mixed-frequency data. Instead, they applied their methods to pre-selected low-dimensional sets of predictors and lags, bypassing the complexity of high-dimensional data structures.

Let us now address the challenge of parameter proliferation, which arises from both the high-dimensionality and the mixed-frequency nature of the data. Several studies have used LASSO as a selection tool to predict corporate bankruptcy; see, for example, petropoulos2020predicting, barbaglia2023forecasting. However, these studies do not address censoring, lack theoretical results, and do not utilize the MIDAS framework. The application of MIDAS in a logistic regression framework for corporate bankruptcy prediction was explored by audrino2019predicting. While their work incorporates the MIDAS approach, it is limited to a low-dimensional set of predictors and does not consider censoring. More closely related to our study, jiang2021predicting examined a penalized logistic regression with the norm $\ell_1$. However, their approach does not account for censoring, lacks theoretical underpinnings, and employs what is referred to as unrestricted MIDAS. Unlike our approach, which uses a restricted MIDAS procedure, unrestricted MIDAS includes all lags as predictors, effectively bypassing the dimension reduction benefits of the MIDAS framework. \footnote{It is worth noting that the term “unrestricted MIDAS" is somewhat misleading, as this approach directly incorporates all lags as independent variables. Consequently, it does not take advantage of the dimension-reduction capabilities inherent in the MIDAS methodology.}

Finally, we compare our theoretical results to the existing literature \textcolor{black}{for both the penalized estimator and its de-sparsified version. The theory of penalized estimators of the high-dimensional logistic regression model has been extensively studied under various situations. As already mentioned, no existing study allows for censoring or approximation error. meier2008group, 10.1214/009053607000000929, buhlmann2011statistics, van_de_geer_estimation_2016 analyzed the logistic regression model using fixed design or isotropic conditions of sub-Gaussian covariates, which are often unsuitable for financial data. van_de_geer_debias developed a de-sparsified inferential procedure for logistic regression with LASSO under uniformly bounded covariates. More recently, caner2023generalized relaxed these assumptions, allowing for random covariate designs with non-normal covariates in the context of penalized Generalized Linear Models (GLM), and established asymptotic theory for the associated de-sparsified estimator. However, compared to the present paper, this work imposes additional assumptions on the shape of the second-order partial derivatives of the loss function. Similarly, han2023high developed the theory for GLM with LASSO by establishing local restricted strong convexity of the loss function, which is related to the quadratic margin condition in the present paper; see Online Appendix (ref). However, they did not propose an inferential procedure.} In the context of mixed-frequency data, babii2022machine, babii2023machine developed the theoretical foundation for high-dimensional time series and panel data linear regression models while accounting for the MIDAS approximation error. However, the theory for penalized logistic regression models incorporating such approximation errors remains unexplored, and none of the aforementioned studies have addressed the challenges posed by censored data.

\paragraph{Outline.} The paper is organized as follows. In Section (ref), we first present the model, followed by a discussion on employing the MIDAS weighting technique and incorporating group structure information among variables. \textcolor{black}{Section (ref) analyzes the estimation properties of the proposed estimator, introduces the associated inferential procedure, and establishes the asymptotic properties of the de-sparsified estimator. Section (ref) presents simulation studies for both prediction and inference.} \textcolor{black}{In Section (ref), we build a dataset on Chinese firm distress to assess the predictive performance of the proposed method, compare it with other approaches, and test each financial covariate for significance.}

\paragraph{Notation.} For $\ell \in \mathbb{N}$, we define $[\ell]=\{1,2, \ldots, \ell\}$. For a vector $\boldsymbol{b} \in \mathbb{R}^p$, its $\ell_q$ norm is denoted as $|\boldsymbol{b}|_q=\left(\sum_{j \in[p]}\left|b_j\right|^q\right)^{1 / q}$ if $q \in[1, \infty)$, \textcolor{black}{its $\ell_{\infty}$-norm is denoted $|\boldsymbol{b}|_{\infty}=\max _{j \in[p]}\left|b_j\right|$, and its empirical $\ell_2$-norm is denoted $\|\boldsymbol{b}\|_N = \sqrt{\left|\boldsymbol{b}\right|_2^2/N}$. Let $diag(b_1, \ldots, b_p)$ be the $p\times p$ diagonal matrix, with diagonal entries equal to $b_1, \ldots, b_p$.} For a matrix $\boldsymbol{A}$, let $\boldsymbol{A}^{\top}$ be its transpose and $\lambda_{\min }(\boldsymbol{A})$ be its smallest eigenvalue. The cardinality of a set $S$ is $|S|$, and we use $S^c$ as its complement. For a vector $\Delta \in \mathbb{R}^p$ and a subset $J \subset[p]$, let $\Delta_J$ be a vector in $\mathbb{R}^p$ with the same coordinates as $\Delta$ on $J$ and zero coordinates on $J^c$, where $J^c$ is the complement of the subset $J$. \textcolor{black}{Let $\Delta_{(J)} := (\Delta_j)_{j \in J} \in \mathbb{R}^{|J|}$ denote the subvector of $\Delta$ corresponding to indices in $J$. For a matrix $\boldsymbol{M} \in \mathbb{R}^{p \times p}$, we define $\boldsymbol{M}_J \in \mathbb{R}^{|J| \times p}$, which is the submatrix consisting of the rows of $\boldsymbol{M}$ corresponding to indices in $J$.} For $a, b \in \mathbb{R}$, we put $a \vee b=\max \{a, b\}$ and $a \wedge b=\min \{a, b\}$. We write $a_N \lesssim b_N$ if there exists a (sufficiently large) absolute constant $v$ such that $a_N \leq v b_N$ for all $N \geq 1$ and $a_N \sim b_N$ if $a_N \lesssim b_N$ and $b_N \lesssim a_N$. The indicator function is denoted by $\mathbbm{1}\{\cdot\}$.

High-dimensional censored MIDAS logistic regression

Logistic regression model

In corporate survival analysis, we focus on the survival time $T$ of a firm. The random variable $T$ represents the duration from the firm's Initial Public Offering (IPO) date to the occurrence of financial distress. Specifically, the IPO date refers to the first day the firm’s stock is publicly traded. Since companies are not listed immediately after their creation, the survival time $T$ in our context differs slightly from the typical survival time considered in traditional survival analysis li2023corporate.

Our main objective is to predict the probability that a firm will survive up to $t$ years, given that it has already been publicly listed for $s$ years. In practice, the survival time $T$ is right-censored by the censoring time $C$, which denotes the duration between the IPO date and the censoring event, occurring at the end of the follow-up period.\footnote{\textcolor{black}{In the empirical application, all the randomness of $C$ comes from the fact that the IPO dates are random and the end date of the follow-up period is fixed.}} Hence, we do not directly observe $T$, but rather the censored value $\widetilde{T} = T \wedge C$, with the indicator $\delta = \mathbbm{1}\{T \leq C\}$.

The financial distress status, indicated by $\mathbbm{1}\{T \leq t\}$, is influenced by covariates $\boldsymbol{Z} \in \mathbb{R}^{K_z}$. We assume, for the moment, that $\boldsymbol{Z}$ has finite variance to ensure well-defined expectations and model the survival indicator $\mathbbm{1}\{T \leq t\}$ using a logistic regression model:

where $\boldsymbol{Z}$ is the covariate vector, and $\boldsymbol{\theta}_0(t,s) \in \mathbb{R}^{K_z}$ is the vector of true parameters specific to $t$ and $s$. For convenience, we use $\boldsymbol{\theta}_0$ as shorthand for $\boldsymbol{\theta}_0(t,s)$. Note that in practice, the model includes an intercept term, which enters the variable $\boldsymbol{Z}$.

Model (ref) is typically estimated via maximum likelihood estimation. This method is based on the characterization of $\boldsymbol{\theta}_0$ as the solution to the population conditional maximum likelihood problem:

In (ref), we have rewritten the classical logistic model's likelihood in a simplified form; see Lemma (ref) for a proof.\footnote{The characterization (ref) is valid under a full-rank condition stated in Lemma (ref).}

However, equation (ref) cannot be directly used for estimation due to the fact that the survival time $T$ is not always observed. To address this issue, we apply the outcome-weighted inverse probability of censoring weighting (OIPCW) method, as outlined by blanche2023logistic.\footnote{An alternative approach for addressing censoring is Inverse Probability Weighting (IPW) bbbe4e97-7833-3faa-820c-4c61f82fe965, zheng2006application, beyhum2024instrumental, beyhum2024dynamic. Further details on both OIPCW and IPW can be found in blanche2023logistic.} This method relies on two standard assumptions about the censoring mechanism, which we describe below. The first assumption is the assumption of independent censoring:

This is a standard assumption in survival analysis. We argue that Assumption (ref) is reasonable in corporate survival analysis because the censoring time for a firm is solely determined by the observation period, with no firms censored before. The second assumption concerns sufficient follow-up:

This assumption implies that some firms have been observed for more than $t$ years without experiencing financial distress, which is necessary for model estimation.

Under Assumptions (ref) and (ref), we obtain an alternative characterization of $\boldsymbol{\theta}_0$, relying only on observed or estimable quantities:

where $H(u) = P(C \geq u|C\ge s)$ is the survival probability of $C$ at time $u$ conditional on $C\ge s$ and $\delta(t) = \mathbbm{1}\{C \geq t \wedge T\}$ is the observation indicator. Equation (ref) is proven in Online Appendix (ref), see Lemma (ref). Essentially, the expectation in (ref) weighs the uncensored observations that fail between $s$ and $t$ by the weights $1/H(t \wedge \widetilde{T})$ to ensure they are representative of firms with survival times between $s$ and $t$.\footnote{The expectation in (ref) is well-defined since $H(t \wedge \widetilde{T})\ge H(t) = P(C \geq t \mid C\ge s) = P\left(C \geq t\right)/P(C\ge s) \ge P(\widetilde{T} \geq t)/P(C\ge s) > 0$ almost surely by Assumption (ref).}

The function $H$ is not directly observed but can be estimated under Assumption (ref) using the classical Kaplan-Meier estimator kaplan1958nonparametric, as described below.

Estimation with mixed-frequency data

Consider an i.i.d. sample of firms \(\left(\widetilde{T}_i, \delta_i, \boldsymbol{Z}_i\right), i \in [N]\), such that for all \(i \in [N]\), \(\widetilde{T}_i \geq s\), i.e., all firms in the sample are observed for at least \(s\) years.\footnote{As shown in the previous section, considering firms with at least \(s\) years of observation does not introduce selection bias under the independent censoring assumption.}

For prediction, suppose we have $K$ covariates and for each covariate $k\in [K]$ we have \(s\) years of lagged covariates \(\left\{z_{i, s - \frac{j-1}{m}, k}: j \in [d]\right\}\), where \(d = s \times m\) represents the total number of lags of each covariate and \(m\) is the annual sampling frequency of the covariate observations. The covariates can be observed at varying frequencies, and although not all lags may enter the regression, we omit such cases for simplicity. The \(k\)-th covariate and its lags are represented as $$ \textcolor{black}{\widetilde{\boldsymbol{Z}}_{i,k} = \left(z_{i, s, k}, z_{i, s - \frac{1}{m}, k}, \ldots, z_{i, s - \frac{d-1}{m}, k}\right)^{\top} \in \mathbb{R}^{d}, k \in [K].} $$ Then the complete vector of lagged covariates and intercept is denoted as $$ \boldsymbol{Z}_i = \left(1,\widetilde{\boldsymbol{Z}}_{i,1}^\top, \widetilde{\boldsymbol{Z}}_{i,2}^\top, \dots, \widetilde{\boldsymbol{Z}}_{i,K}^\top\right)^\top\in\mathbb{R}^{K_z} $$ where \(K_z = K \times d + 1\). The function \(H\) in (ref) can be estimated using the Kaplan-Meier estimator: $$ \widehat{H}(u) = \prod_{j \leq u} \left(1 - \frac{dN(j)}{\widetilde{T}(j)}\right), $$ where \(N(j) = \sum_{i=1}^N \mathbbm{1}\{\widetilde{T}_i \leq j, \delta_i = 0\}\) is the counting process for censoring events up to time $j$, and \(dN(j) = N(j) - \lim_{j' \to j, j' < j} N(j')\) denotes the jump of the process \(N\) at time \(j\). Additionally, \(\widetilde{T}(j) = \sum_{i=1}^N \mathbbm{1}\{\widetilde{T}_i \geq j\}\) is the number of units who are under observation just before time \(j\).

We consider datasets that are high-dimensional. For instance, in our empirical application, as summarized in Table (ref), if we consider firms that have survived \(s = 6\) years, with \(K = 95\) covariates measured \(m = 4\) times per year, the total number of parameters to estimate is \(6 \times 4 \times 95 + 1 = 2,281\), including the intercept. When the sample size is not much larger than the number of parameters, the curse of dimensionality arises, complicating computations and reducing estimation precision.

To address this, dimension-reduction techniques are necessary. A common approach is to directly apply the LASSO tibshirani1996regression to the original predictors. For an i.i.d. sample \(\{(\widetilde{T}_i, \delta_i, \boldsymbol{Z}_i), i \in [N]\}\), the $\ell_1$-norm penalized estimator solves the empirical version of (ref) with the LASSO penalty.

Here, we follow a different approach babii2022machine to reduce the dimension based on Mixed-Data Sampling ghysels2006predicting, which is designed to address parameter proliferation in mixed-frequency data. MIDAS approximates the coefficients of high-frequency lag polynomials using a finite dictionary of functions. Specifically, let us write

where \(\omega_k:[0,1]\mapsto \mathbb{R}\), \(k \in [K]\), are weight functions for the lag polynomials such that $ \omega_k\left(\frac{j-1}{d}\right)/d=\theta_{0,1+d (k-1) +j}$. Let \(\{w_l: l = 1, \ldots, L\}\) be the dictionary of functions \textcolor{black}{and $L$ is the size of the dictionary.} For each \(k \in [K]\), we assume there exist coefficients \(\boldsymbol{\beta}_{0,k}^* = (\beta_{0,k,1}^*, \beta_{0,k,2}^*,\ldots, \beta_{0,k,L}^*)^\top \in \mathbb{R}^L\) such that: $$ \omega_k(u) \approx \sum_{l=1}^L \beta_{0,k,l}^* w_l(u), \quad u \in [0, 1]. $$ This reduces the number of parameters from \(K \times d+1 \) to \(K \times L+1 \). The simplest dictionary consists of algebraic power polynomials (e.g., Almon polynomials almon1965distributed), but other orthogonal bases of \(L_2[0,1]\) can be used to improve performance with correlated covariates.\footnote{\(L_2[0,1]\) denotes the space of square-integrable functions \(f: [0,1] \to \mathbb{R}\).} {\color{black}With the MIDAS approximation, we define the following MIDAS-weighted covariate vector \[ \boldsymbol{X}_i = \left(1,\widetilde{\boldsymbol{Z}}_{i,1}^\top W, \widetilde{\boldsymbol{Z}}_{i,2}^\top W, \ldots, \widetilde{\boldsymbol{Z}}_{i,K}^\top W\right)^\top \in \mathbb{R}^{p}, \] where $p = K \times L + 1$ is the dimension, \(W = \left(w_l\bigl(\frac{j-1}{d}\bigr)/d\right)_{j \in [d], l \in [L]} \in \mathbb{R}^{d \times L}\) is the weighting matrix, and we refer to $W^{\top}\widetilde{\boldsymbol{Z}}_{i,k}$ as the $k$-th group of MIDAS-weighted covariates. Although the number of parameters to be estimated is reduced by taking a small MIDAS dictionary size $L$, a large $K$ can still pose a high-dimensional challenge. }

{\color{black} To further reduce dimensionality and exploit the natural grouping of covariates, we employ the sparse-group LASSO simon2013sparse. While standard LASSO induces sparsity at the individual-variable level and group LASSO enforces sparsity at the group level, sparse-group LASSO provides a more flexible framework by encouraging sparsity both within and across groups. This is achieved through a convex combination of the LASSO and group LASSO penalties. As a result, sparse-group LASSO can select relevant groups while simultaneously performing variable selection within those groups, making it particularly well suited to high-dimensional settings with hierarchical covariate structures, such as in high-dimensional MIDAS regressions. For \(\boldsymbol{\beta}\in\mathbb{R}^{p}\), define the sparse-group LASSO penalty as $$ \Omega(\boldsymbol{\beta}) = \alpha |\boldsymbol{\beta}|_1 + (1-\alpha) \|\boldsymbol{\beta}\|_{2,1}, \quad \|\boldsymbol{\beta}\|_{2,1} = \sum_{G \in \mathcal{G}} |\boldsymbol{\beta}_G|_{2}, $$ where $0 \leq \alpha \leq 1$ is the weight parameter that balances the LASSO and group LASSO penalties, and $\mathcal{G}$ denotes a pre-specified collection of groups, with each group $G \in \mathcal{G}$ being a subset of $[p]$.} {\color{black} In our analysis, the specified group structure is \(\mathcal{G} = \{G_k : k \in [K+1]\}\), where $G_1=\{1\}$ corresponds to the intercept and for $k = 2, \ldots, K+1$, \(G_k = \{(1+(k-2)L)+ 1, \ldots, 1+(k-1)L\}\) corresponds to the $(k-1)$-th group of MIDAS-weighted covariates $W^{\top}\widetilde{\boldsymbol{Z}}_{i,(k-1)}$. This grouping arises by treating each covariate together with all of its lags as a single group, and then applying the MIDAS weighting matrix $W$ to each such group of lags. As shown by babii2022machine, this simple and natural group structure enhances prediction performance. Alternative group structures are also possible. For example, forming groups by pairing conceptually related covariates, such as Return on Assets (ROA) and Return on Equity (ROE). We explore the robustness of our results to such alternative groupings in the empirical application.}

\textcolor{black}{For ease of reference, the notation introduced so far is summarized in Table (ref).}

Finally, we propose the sparse-group LASSO estimator, which minimizes

where $$ R_N(\boldsymbol{\beta}) = \frac{1}{N} \sum_{i=1}^N -\frac{\delta_i(t) \mathbbm{1}\{\widetilde{T}_i \leq t\}}{\widehat{H}(t \wedge \widetilde{T}_i)} \boldsymbol{X}_i^\top \boldsymbol{\beta} + \log \left(1 + \exp(\boldsymbol{X}_i^\top \boldsymbol{\beta})\right). $$ Here, \(\lambda \geq 0\) controls the regularization, and \(\alpha \in [0,1]\) balances the sparsity and the group structure. Choosing \(\alpha = 1\) recovers LASSO, while \(\alpha = 0\) corresponds to group LASSO.\footnote{In our practical implementation, we do not penalize the intercept coefficient. For simplicity, we do not write this in the equations.} \textcolor{black}{In practice, $\alpha$ can be fixed or selected jointly with $\lambda$ by using cross-validation.} We call sg-LASSO-MIDAS the approach embodied by (ref).

In the simulation and empirical sections, we also examine two alternative methods as benchmarks. The first method LASSO-UMIDAS employs unrestricted lag polynomials combined with LASSO. Unlike the MIDAS approach, LASSO-UMIDAS does not impose restrictions on the polynomials, resulting in the need to estimate a significantly larger number of parameters. Moreover, no structural constraints are applied to the coefficients of the lags. As a second alternative, we consider LASSO-MIDAS, which adopts the MIDAS approximation employed in the sg-LASSO-MIDAS method but avoids the group penalty. It corresponds to (ref) with a fixed mixing parameter of $\alpha = 1$.

Theoretical results

\textcolor{black}{In this section, we first outline the main assumptions for the proposed penalized estimator (ref) and analyze the finite sample properties of the sparse-group LASSO estimator. Then we introduce the de-sparsified procedure and discuss the asymptotic properties of the de-sparsified sparse-group LASSO estimator. Both the (de-sparsified) LASSO and (de-sparsified) group LASSO estimators are covered as special cases.\footnote{We treat $\alpha$ as constant for the theory but optimize it through cross-validation in practice.}} Recall that we have an i.i.d. sample $\{(\widetilde{T}_i,\delta_i,\boldsymbol{Z}_i),\ i\in[N]\}$ such that $\widetilde{T}_i\ge s$ for all $i\in[N].$ We consider an asymptotic regime where $N$ goes to infinity and $p$ goes to infinity as a function of $N$. High-dimensional $\ell_1$-norm penalized logistic regression has been studied in the literature; see, for instance, 10.1214/009053607000000929,van_de_geer_estimation_2016,caner2023generalized and han2023high. However, none of these studies account for censoring or allow for approximation errors. Instead, we explicitly take into account the approximation error stemming from the MIDAS approximation defined as $$ E_i = \boldsymbol{Z}_{i}^{\top}\boldsymbol{\theta}_0 - \boldsymbol{X}_{i}^{\top}\boldsymbol{\beta}_0, \quad i \in [N], $$ where $\boldsymbol{\beta}_0 = \left(\theta_{0,1}, (\boldsymbol{\beta}_{0,1}^*)^{\top},(\boldsymbol{\beta}_{0,2}^*)^{\top},\ldots,(\boldsymbol{\beta}_{0,K}^*)^{\top}\right)^{\top}\in \mathbb{R}^p$ is the true parameter of interest \textcolor{black}{and we assume it is exactly sparse. We refer to $\theta_{0,1}$ as the intercept term and $\boldsymbol{\beta}_{0,k}^* \in \mathbb{R}^{L}$ as the $k$-th group of MIDAS-weighted parameters. Let $\boldsymbol{E} = (E_1, E_2, \ldots, E_N)^{\top}$ collect all approximation errors.}

\textcolor{black}{If one is instead interested in allowing for approximate sparsity of $\boldsymbol{\beta}_{0}$ van_de_geer_estimation_2016,belloni2018high,babii2022machine, a straightforward modification is to replace $E_i$ by the following error term $\widetilde{E}_i$ defined by \[ \widetilde{E}_i = \underbrace{\boldsymbol{Z}_i^{\top}\boldsymbol{\theta}_0 - \boldsymbol{X}_i^{\top}\boldsymbol{\beta}_0}_{\text{MIDAS approximation error}} +\; \underbrace{\boldsymbol{X}_i^{\top}\boldsymbol{\beta}_0 - \boldsymbol{X}_i^{\top}\boldsymbol{\beta}}_{\text{approximate sparsity error}}, \] where $\boldsymbol{\beta}$ is a candidate oracle that is exactly sparse, and the assumptions (ref) and (ref) below should therefore be imposed on $\boldsymbol{\beta}$ rather than on $\boldsymbol{\beta}_0$.}\footnote{\textcolor{black}{Alternatively, by adapting the argument underlying Corollary 12.7 of van_de_geer_estimation_2016, one can also derive a finite-sample bound that yields a separate treatment of the MIDAS approximation error and the approximate sparsity error.}} \textcolor{black}{In this paper, we stick to exact sparse $\boldsymbol{\beta}_0$ to simplify the theoretical analysis.}

Estimation theory

We start by introducing the following assumptions.

This condition just requires that the variables have more than $4$ finite moments, allowing for polynomial tails commonly observed with financial variables. \textcolor{black}{It is worth noting that Assumption (ref) also implicitly controls correlations among covariates, as extremely high correlations would lead some linear combination to have an excessively large $q$-th moment, thereby violating the bound.}

Assumption (ref) is similar to the compatibility condition discussed in 10.1214/009053607000000929,van_de_geer_estimation_2016,caner2023generalized, as well as the restricted Fisher-information matrix eigenvalue condition described in han2023high. This is a high-dimensional version of the full-rank condition guaranteeing the asymptotic properties of the maximum likelihood estimator in the (low-dimensional) logistic regression with misspecification error.

Next, we need to introduce additional definitions. Let $S_{\boldsymbol{\beta}_0}=\{j\in[p]:\ \boldsymbol{\beta}_{0,j}\ne 0\}$ and $ \mathcal{G}_{\boldsymbol{\beta}_0}=\{G \in \mathcal{G}:\ \left(\boldsymbol{\beta}_{0}\right)_{G}\ne \boldsymbol{0}\}$ be the support and the group support of the target parameter $\boldsymbol{\beta}_0$. Let $\sqrt{s_{\boldsymbol{\beta}_0}} = \alpha \sqrt{|S_{\boldsymbol{\beta}_0}|} + (1 - \alpha)\sqrt{|\mathcal{G}_{\boldsymbol{\beta}_0}|}$ be the sparsity level and $G^*=\max _{G \in \mathcal{G}_{\boldsymbol{\beta}_0}}|G|$ be the size of the largest group in $\mathcal{G}_{\boldsymbol{\beta}_0}$. For simplicity, we suppose that $s_{\boldsymbol{\beta}_0}\ge 1$ (otherwise, it suffices to replace $s_{\boldsymbol{\beta}_0}$ by $s_{\boldsymbol{\beta}_0}\vee 1$ in all assumptions and bounds involving $s_{\boldsymbol{\beta}_0} $). We impose the following assumption on $s_{\boldsymbol{\beta}_0} $.

This is a condition on the degree of sparsity $s_{\boldsymbol{\beta}_0} $, the size of the largest group $G^*$, the $\ell_1$-norm of the approximation error $|\boldsymbol{E}|_1 = \sum_{i = 1}^{N} |E_i|$, and the relative growth rate of $N$ and $p$. The condition is more likely to hold when $s_{\boldsymbol{\beta}_0} $, $G^*$, $|\boldsymbol{E}|_1$ or $1/q$ are smaller and $p$ does not grow too quickly with $N$. Note also that if $\lambda$ is too low or too large, the condition might fail to hold. This condition allows establishing a connection between empirical and population effective sparsity, enabling the extension of the quadratic margin condition to its sampled version. 10.1214/09-EJS506 briefly discussed it specifically for data with Gaussian tails in the case of the LASSO for the linear model. We extend this framework to accommodate heavy-tailed data and approximation error in the logistic regression model. When there is no such approximation error, that is $E_i = 0, $ for all $i \in [N]$, a similar assumption imposed on $\lambda$ and $s_{\boldsymbol{\beta}_0}$ is used in van_de_geer_estimation_2016 and han2023high.

We now establish bounds on the estimation error, presenting two distinct types. The first type pertains to the parameter estimation error $\Omega\left(\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta}_0\right)$, while the second focuses on prediction accuracy. Consider a scenario where, for some $\boldsymbol{z} = \left(1, \tilde{\boldsymbol{z}}_1, \ldots, \tilde{\boldsymbol{z}}_K \right)^{\top} \in \mathbb{R}^{K_z}$, where $\tilde{\boldsymbol{z}}_k \in \mathbb{R}^{d}, k \in [K]$, we aim to estimate $P(\boldsymbol{z}) = P(T \leq t \mid \boldsymbol{Z} = \boldsymbol{z}, T \geq s)$, representing the probability that a firm with covariates $\boldsymbol{z}$, having survived at least $s$ years, fails before $t$. We estimate $P(\boldsymbol{z})$ using $\widehat{P}(\boldsymbol{z}) = \frac{\exp \left( \boldsymbol{x}^{\top} \widehat{\boldsymbol{\beta}} \right)}{1 + \exp \left(\boldsymbol{x}^{\top}\widehat{\boldsymbol{\beta}}\right)}$ where $\boldsymbol{x} = \left(1, \tilde{\boldsymbol{z}}_1^{\top}W, \ldots, \tilde{\boldsymbol{z}}_K^{\top}W\right)^{\top}$. Our goal is to provide a bound for the error $\widehat{P}(\boldsymbol{z}) - P(\boldsymbol{z})$. This bound will depend on the term $e = \boldsymbol{z}^\top \boldsymbol{\theta}_0 - \boldsymbol{x}^\top \boldsymbol{\beta}_0$, which represents the MIDAS approximation error at the covariate $\boldsymbol{z}$. The following theorem formally states this result.

Let us now discuss the theorem. First, we require that $\lambda$ at least has the same order with respect to $ p^{\frac{1}{q}} \sqrt{\log p}/N^{\frac{1}{2} - \frac{1}{q}}. $ For bounds on the LASSO under sub-Gaussian errors, it suffices that \(\lambda\) is of the order of \(\sqrt{\log p / N}\). Our condition is stricter due to the presence of heavy-tailed variables. However, as \( q \to \infty \), the variables are no longer heavy-tailed, and we recover the order \(\sqrt{\log p / N}\) for \(\lambda\). As is standard in the literature, in practice, we select \(\lambda\) in practice via cross-validation (see Sections (ref) and (ref)). The dependence of our rates on \(\lambda\) aligns with those for the standard LASSO estimator in high-dimensional regression.

Our bounds also depend on the \(\ell_1\)-norm \(|\boldsymbol{E}|_1\) of the approximation error. To the best of our knowledge, this work is the first to establish results for the high-dimensional logistic regression model with an approximation error. Such results, however, are well-established for the LASSO in linear regression bickel2009simultaneous. To achieve this result, we bound the difference between the empirical and population loss functions not only by terms related to the empirical process but also by a term dependent on the approximation error \(\boldsymbol{E}\). Addressing this challenge is particularly difficult due to the nonlinearity of the problem. Interested readers are referred to Online Appendix (ref) for the detailed proof of Theorem (ref). Regarding the estimated prediction probability, its error bound matches the order of the parameter estimation error, with an additional term which is a function of the MIDAS approximation error \( e \).

As a concluding remark, we note that the Kaplan–Meier estimator converges at the $1/\sqrt{N}$ rate. Hence, estimating the weights $1/H(t \wedge \widetilde{T}_i)$ does not affect the convergence rate of our penalized estimator.

{\color{black} We then consider a regime where the MIDAS approximation error vanishes at a specific rate, as described in the following assumption.

The following corollary on parameter estimation follows immediately.

Corollary (ref) characterizes the impact of censoring and heavy-tailed covariates on the estimation accuracy of $\widehat{\boldsymbol{\beta}}$. Unlike Corollary 3.1 of babii2022machine, which explicitly incorporates persistence in the dependent variable for a linear time-series model, our setting involves i.i.d. censored data in a logistic regression framework with lagged covariates treated as separate variables. In this context, the $q$-th moment bound in Assumption (ref) implicitly controls the effect of correlations among lagged covariates. Hence, Corollary (ref) captures the impact of covariate correlation strength indirectly. For literature that explicitly incorporates covariate correlations in error bounds for i.i.d. data, we refer to hebiri2012correlations,lassoprediction. While these works focus on linear models, their arguments could potentially be extended to nonlinear settings, which we leave for future research. } {\color{black}

De-sparsified procedure and inference theory

Beyond estimation with the proposed sparse-group LASSO estimator, our second goal is to conduct statistical inference for the underlying parameters. Let us use the notation $ \boldsymbol{\beta}_0 := \left(\beta_{0,1}, \beta_{0,2}, \ldots, \beta_{0,p}\right)^{\top} \in \mathbb{R}^p. $ We focus on inference on the subvector $$ \boldsymbol{\beta}_{0,(\mathcal{J})} := (\beta_{0,j})_{j \in \mathcal{J}} \in \mathbb{R}^{|\mathcal{J}|}, $$ where $\mathcal{J} \subseteq [p]$ denotes a fixed index set. As an example, one may wish to test the hypothesis $\mathbb{H}_0: \boldsymbol{\beta}_{0,(\mathcal{J})} = \boldsymbol{0}$. However, it is well known that the penalized estimator $\widehat{\boldsymbol{\beta}}$ is asymptotically biased. Hence, a debiasing (or desparsification) procedure is required.

Our de-sparsified methodology builds on the approach of van_de_geer_debias, which approximates the inverse of the population Gram matrix $\boldsymbol{\Sigma}_{\boldsymbol{\beta}_0}$ via nodewise regressions and employs this approximation to correct the bias of the penalized estimator, yielding de-sparsified estimates with asymptotically normal components. Unlike in the linear model, the empirical Gram matrix in the logistic model depends on the penalized estimator $\widehat{\boldsymbol{\beta}}$. In our setting, the sparse-group LASSO estimator $\widehat{\boldsymbol{\beta}}$ in (ref) is further complicated by censoring, MIDAS approximation errors, and the use of a sparse-group LASSO penalty, distinguishing our approach from the existing literature.

We now describe the procedure to estimate $\boldsymbol{\Sigma}_{\boldsymbol{\beta}_0}^{-1}$. To this end, we define the population Gram matrix in the presence of the MIDAS approximation error $\boldsymbol{E}$:

Let $\boldsymbol{\Theta} = (\boldsymbol{\Theta}_1,\ldots,\boldsymbol{\Theta}_p)^{\top} = \boldsymbol{\Sigma}_{\boldsymbol{\beta}_0}^{-1}$, and let $\widehat{\boldsymbol{\Theta}} = \left(\widehat{\boldsymbol{\Theta}}_1, \ldots, \widehat{\boldsymbol{\Theta}}_p\right)^{\top}$ denote its estimator. Each $\widehat{\boldsymbol{\Theta}}_j \in \mathbb{R}^p$, $j \in [p]$, can be obtained by using a weighted nodewise regression. The weighting arises naturally from the logistic model, as the Gram matrix in logistic regression incorporates observation-specific weights. Accordingly, we define the weights used in the nodewise regressions at both the population and sample levels. $$

$$ where both are $N \times N$ diagonal matrices, with

For $\boldsymbol{\beta} \in \mathbb{R}^{p}$, which may represent either $\boldsymbol{\beta}_0$ or $\widehat{\boldsymbol{\beta}}$, let $\boldsymbol{X}_{\boldsymbol{\beta}}:=W_{\boldsymbol{\beta}} \boldsymbol{X} \in \mathbb{R}^{N \times p}$, and denote the $j$-th column of that matrix as $\boldsymbol{X}_{\boldsymbol{\beta}, j} \in \mathbb{R}^{N}$, and the matrix containing all the columns of $\boldsymbol{X}_{\boldsymbol{\beta}}$ except the $j$-th one as $\boldsymbol{X}_{\boldsymbol{\beta},-j} \in \mathbb{R}^{N \times (p-1)}$.

To implement the weighted nodewise regressions, for each $j \in [p]$ we solve

where $ \widehat{\boldsymbol{\gamma}}_{\widehat{\boldsymbol{\beta}}, j} = ( \widehat{\gamma}_{\widehat{\boldsymbol{\beta}}, j,1}, \ldots, \widehat{\gamma}_{\widehat{\boldsymbol{\beta}}, j,j-1},\widehat{\gamma}_{\widehat{\boldsymbol{\beta}}, j,j+1},\ldots, \widehat{\gamma}_{\widehat{\boldsymbol{\beta}}, j,p})^{\top} \in \mathbb{R}^{p-1}$ and $\lambda_j>0$ is the penalty parameter for each regression, which can be selected using cross-validation in practice. Then the estimator of $\boldsymbol{\Theta}$ is $$ \widehat{\boldsymbol{\Theta}} = \left(\widehat{\boldsymbol{\Theta}}_1, \ldots, \widehat{\boldsymbol{\Theta}}_p\right)^{\top} = \widehat{\boldsymbol{B}}^{-1}\widehat{\boldsymbol{C}}, $$ where

Each $\widehat{\boldsymbol{\Theta}}_j$ can also be written as

To establish the convergence rate of $\widehat{\boldsymbol{\Theta}}$, we introduce the following notation and assumptions. For each $j \in [p]$, consider the population nodewise regressions

where $ \boldsymbol{\gamma}_{\boldsymbol{\beta}_0, j} = ( \boldsymbol{\gamma}_{\boldsymbol{\beta}_0, j, 1}, \ldots,\boldsymbol{\gamma}_{\boldsymbol{\beta}_0, j, j-1},\boldsymbol{\gamma}_{\boldsymbol{\beta}_0, j, j+1},\ldots \boldsymbol{\gamma}_{\boldsymbol{\beta}_0, j, p-1})^{\top} \in \mathbb{R}^{p-1}$ and define the error as

Assumption (ref) (i) requires that the MIDAS approximation errors vanish sufficiently fast asymptotically.\footnote{\textcolor{black}{If each row of $\boldsymbol{\Theta}$ has a bounded $\ell_1$-norm, we only require $\|\boldsymbol{E}\|_{N}=o_{p}(N^{-1/2})$ in Assumption (ref) (i).}} Assumption (ref) (ii) imposes a stronger moment condition on the covariates than Assumption (ref), requiring $q>12$. This requirement differs from much of the existing literature. For example, Condition (D4) in van_de_geer_debias and Assumption 3 in caner2023generalized either assume uniformly bounded covariates or impose a uniform lower bound on $w^2_{\boldsymbol{\beta}_0,i} \ge c$ for all $i \in [N]$, where $c$ is a strictly positive constant. Such assumptions may be restrictive or unrealistic when covariates have heavy tails. In contrast, our framework allows for heavy-tailed covariates without requiring boundedness or uniform positivity of $w^2_{\boldsymbol{\beta}_0,i}$. This added generality comes at the cost of a more restrictive sparsity condition, as formalized in Assumption (ref). Assumption (ref) (iii) can be interpreted as a weaker version of the second statement of Condition (D1) in van_de_geer_debias, where the latter needed $|\boldsymbol{\eta}_{\boldsymbol{\beta}_0,j}|_{\infty}$ to be bounded. Assumption (ref) (iv) weakens Assumption 4.1 (ii) in babii2023machine, which imposes that every row of $\boldsymbol{\Theta}$ has a uniformly bounded $\ell_1$-norm. Together with Assumption (ref) (ii), Assumption 4.1 (ii) in babii2023machine implies our Assumption (ref) (iv), so our requirement is less restrictive.\footnote{\textcolor{black}{If each row of $\boldsymbol{\Theta}$ has a bounded $\ell_1$-norm, we have $\mathbbm{E}(|\boldsymbol{\Theta}_{j}^{\top}\boldsymbol{X}_i|^{q}) = |\boldsymbol{\Theta}_{j}|^q_2\mathbbm{E}(|\boldsymbol{\Theta}_{j}^{\top}\boldsymbol{X}_i/|\boldsymbol{\Theta}_{j}|_2|^{q}) \leq |\boldsymbol{\Theta}_{j}|^q_1K_0$ where the inequality follows from Assumption (ref) (ii) and the fact that $|\boldsymbol{\Theta}_{j}|_2 \leq |\boldsymbol{\Theta}_{j}|_1$.}}

Define $\bar{s} = \max_{j \in \mathcal{J}} |S_j|$, where $S_j$ denotes the index set of nonzero elements of $\boldsymbol{\Theta}_j$, and set $s^{*} = \bar{s} \vee s_{\boldsymbol{\beta}_0}$. We have the following assumption for the sparsity level $s^*$.

When the covariates are bounded, and there is no MIDAS approximation error, the admissible sparsity level for valid inference is typically of the order $o\big(\sqrt{N}/\log p\big)$ van_de_geer_debias. The sparsity condition imposed in this paper is more restrictive, reflecting the presence of MIDAS approximation error and heavy-tailed covariates.

Although lemmas of this type for nodewise regression are not new in the literature, this paper makes several distinct contributions. van_de_geer_debias established a related result for nodewise regression with the LASSO. However, their analysis relies on strong assumptions, including uniformly bounded covariates and uniformly bounded products of nodewise regression coefficients and covariates. In addition, their framework does not accommodate censored data or the MIDAS approximation error. Similarly, caner2023generalized derived a general result for weighted nodewise regression using a structured sparsity estimator. Nevertheless, their analysis neither considers censored data nor accounts for the MIDAS approximation error, and it assumes that $w^2_{\boldsymbol{\beta}_0,i}$ is uniformly bounded away from zero. Such an assumption can be restrictive and is generally violated in the presence of unbounded covariates.

Given the estimator $\widehat{\boldsymbol{\Theta}}$, the de-sparsified sparse-group LASSO estimator $\widehat{\boldsymbol{b}} = (\widehat{b}_1, \ldots, \widehat{b}_p)^{\top}$ is constructed componentwise as

where $\widehat{\beta}_j$ is the $j$-th coordinate of the sparse-group LASSO estimator $\widehat{\boldsymbol{\beta}}$. We call de-sparsified sg-LASSO-MIDAS the approach embodied by (ref). To derive the limiting distribution of the de-sparsified estimator, one more assumption is required.

Assumption (ref) is a typical assumption for high-dimensional logistic regression, see a similar condition (vi) in Theorem $3.3$ of van_de_geer_debias or Assumption $8$ of caner2023generalized, which is used to show Lyapunov’s condition for the limiting distribution of the de-sparsified sparse-group LASSO estimator. We extend it to the case of censored data.

The following theorem establishes the asymptotic distribution of the proposed de-sparsified estimator, which is applicable to high-dimensional censored data. The detailed proof is provided in Online Appendix (ref). In order to state the theorem, we introduce $$ \xi_{KM, i} = \mathbbm{E}\left(\boldsymbol{X}_{(2)}\frac{\delta_{(2)}(t) \mathbbm{1}{\{\widetilde{T}_{(2)} \leq t\}}}{H^2\left(t \wedge \widetilde{T}_{(2)}\right)}v\left(\widetilde{T}_i, \delta_i,t \wedge \widetilde{T}_{(2)}\right)\bigg| \boldsymbol{X}_i, \widetilde{T}_i, \delta_i\right), $$ where $\left(\widetilde{T}_{(2)}, \delta_{(2)}, \boldsymbol{X}_{(2)}\right)$ denotes an independent copy drawn from the same distribution as the sample $\left\{\left(\widetilde T_i, \delta_{i}, \boldsymbol{X}_i\right), i \in [N]\right\}$. Here, $v\left(\widetilde{T}_i, \delta_i, \cdot \right)$ is the influence function of the Kaplan-Meier estimator of $H$ for observation $\left(\widetilde T_i, \delta_i\right)$ lo1986product,gill2006lectures. Specifically, for any positive value $z$, we have

where $\bar{P}(u)=P\left(\widetilde{T}>u\right)$ and $P_0(du)=\frac{\partial P\left(\widetilde{T} \leq u, \delta = 0\right)}{\partial u}du$.

Theorem (ref) establishes the limiting distribution of the de-sparsified sparse-group LASSO estimator $\widehat{\boldsymbol{b}}_{(\mathcal{J})}$. In contrast to traditional de-sparsified estimators studied in van_de_geer_debias, caner2023generalized, and babii2023machine, the variance of $\widehat{\boldsymbol{b}}_{(\mathcal{J})}$ is influenced by censored data, a phenomenon that, to the best of our knowledge, has not been previously studied.

Finally, we provide a plug-in estimator of the variance $\boldsymbol{V}_\mathcal{J}$: \[ \widehat{\boldsymbol{V}}_\mathcal{J} = \widehat{\boldsymbol{\Theta}}_\mathcal{J}\left(\frac{1}{N}\sum_{i=1}^{N}\widehat{\boldsymbol{\sigma}}_i\widehat{\boldsymbol{\sigma}}_i^{\top}\right)\widehat{\boldsymbol{\Theta}}_\mathcal{J}^{\top}, \] where \[ \widehat{\boldsymbol{\sigma}}_i = \left(-\frac{\delta_i(t) \mathbbm{1}\{\widetilde{T}_i \leq t\}}{\widehat{H}(t \wedge \widetilde{T}_i)}+\frac{\exp(\boldsymbol{X}_{i}^{\top}\widehat{\boldsymbol{\beta}})}{1+\exp(\boldsymbol{X}_{i}^{\top}\widehat{\boldsymbol{\beta}}) }\right)\boldsymbol{X}_{i} - \frac{1}{N}\sum_{k=1}^{N}\boldsymbol{X}_{k}\frac{\delta_k(t) \mathbbm{1}{\{\widetilde{T}_k \leq t\}}}{\widehat{H}^2\left(t \wedge \widetilde{T}_k\right)}\hat{v}\left(\widetilde{T}_i, \delta_i,t \wedge \widetilde{T}_k\right), \] and $\hat{v}\left(\widetilde{T}_i, \delta_i, \cdot\right)$ is the plug-in estimator of $v\left(\widetilde{T}_i, \delta_i, \cdot\right)$. Specifically, for any positive value $z$,

where $\bar{P}_N(u) = \frac{1}{N}\sum_{l=1}^N\mathbbm{1}\{\widetilde{T}_l > u\}$. \footnote{\textcolor{black}{For $\hat{v}$, its implementation can be found in the R package 'survival'; see the argument 'influence' of the function 'survfit' in therneau2015package.}} The term $\frac{1}{N}\sum_{j=1}^N\frac{\mathbbm{1}\{\widetilde{T}_j \leq z \wedge \widetilde{T}_i, \delta_j = 0\}}{\bar{P}_N^2(\widetilde{T}_j)}$ in (ref) is an estimate of the term $\int_0^{z \wedge \widetilde{T}_i} \frac{P_0(du)}{\bar{P}^2(u)}$ in (ref) because we have $ \int_0^{z \wedge \widetilde{T}_i} \frac{1}{\bar{P}^2(u)}P_0(du) = \mathbbm{E}\left(\frac{1}{\bar{P}^2(\widetilde{T})}\mathbbm{1}\{\widetilde{T} \leq z \wedge \widetilde{T}_i\}\mathbbm{1}\{\delta=0\}\right), $ since $P_0(du)$ is the probability mass of $\widetilde{T}$ in an infinitesimal interval around $u$, computed only over the event $\delta=0$. }

Simulations

We first evaluate the predictive performance of three methods through simulations: i) LASSO-UMIDAS, an unstructured LASSO estimator \textcolor{black}{with unrestricted MIDAS weights}, ii) LASSO-MIDAS, an unstructured LASSO estimator using MIDAS weights, and iii) sg-LASSO-MIDAS, a structured sparse-group LASSO estimator with MIDAS weights. The sg-LASSO-MIDAS approach, specifically, highlights the advantages of leveraging group structures and dictionaries within a high-dimensional framework, offering a compelling comparison to LASSO-MIDAS and LASSO-UMIDAS babii2022machine. The prediction simulation results showcase the method's strengths in achieving superior prediction accuracy with finite sample data. \textcolor{black}{In addition to predictive performance, we evaluate the inference performance of de-sparsified sg-LASSO-MIDAS. The inference simulation results indicate that it provides accurate finite-sample inference, delivering empirical rejection rates close to the nominal level across various scenarios.}

Simulation design

Let us describe the data-generating process. There are $K=50$ high-frequency covariates, but only the first two enter the model. All the observations in the simulated dataset have survived at least $s$ years, and we are interested in a yearly/quarterly frequency $ m = 4$. We consider $s=6$ years of lagged data.

For the generation of $z_{i, \frac{j}{m}, k}, j \in [d]$, we first initiate the processes by letting $\left(z_{i,\frac{1}{m},1}, \ldots, z_{i,\frac{1}{m},K}\right)^{\top}$ \textcolor{black}{follow a $\mathcal{N}\left(\boldsymbol{0}, \Sigma \right)$} distribution with $\Sigma_{u, v}=\rho_{0}^{|u-v|}, u,v \in [K]$. Then, the high-frequency covariates $z_{i, \frac{j}{m},k}, j \in [d], k \in [K]$ are generated according to the following scenarios:

\newcounter{scenario} \refstepcounter{scenario} \noindentScenario \thescenario: $z_{i, \frac{j}{m},k}=\rho z_{i,\frac{j-1}{m},k}+\nu_{i,k}$, $k \in [K], j \in \{2,3,\ldots,d\}$, and $(\nu_{i,1}, \ldots, \nu_{i,K})^{\top} \sim_{\text {i.i.d }} \mathcal{N}\left(\boldsymbol{0}, \Sigma(1-\rho^2)\right)$ with $\rho=0.1$ and $\rho_0 = 0.1$.

\refstepcounter{scenario} \noindentScenario \thescenario: $z_{i, \frac{j}{m},k}=\rho z_{i,\frac{j-1}{m},k}+\nu_{i,k}$, $k \in [K], j \in \{2,3,\ldots,d\}$, and $(\nu_{i,1}, \ldots, \nu_{i,K})^{\top} \sim_{\text {i.i.d }} \mathcal{N}\left(\boldsymbol{0}, \Sigma(1-\rho^2)\right)$ with $\rho=0.9$ and $\rho_0 = 0.1$.

It is clear that $\rho$ regulates the degree of series dependence among original lagged covariates, while $\rho_0$ represents the level of cross-covariate dependence across all $K$ covariates. \textcolor{black}{ We refer to $\rho$ as the autocorrelation strength of original lags and $\rho_0$ as cross-covariate correlation strength. Inspired by the empirical application, we consider two more scenarios that allow the covariates to have heavy tails. In these cases, we first initiate the processes with $\left(z_{i,\frac{1}{m},1}, \ldots, z_{i,\frac{1}{m},K}\right)^{\top} \sim$ student-$t(2)$ with the covariance matrix $\Sigma_{u, v}=\rho_{0}^{|u-v|}, u,v \in [K]$. Then, the third and fourth scenarios are as follows.}

\refstepcounter{scenario} \noindentScenario \thescenario: $z_{i, \frac{j}{m},k}=\rho z_{i,\frac{j-1}{m},k}+\nu_{i,k}$, $k \in [K], j \in \{2,3,\ldots,d\}$, and $(\nu_{i,1}, \ldots, \nu_{i,K})^{\top} \sim_{\text {i.i.d}}$ student-$t$ with degree $2$ and its covariance matrix $\Sigma(1-\rho^2)$, with $\rho=0.1$ and $\rho_0 = 0.1$.

{\color{black} \refstepcounter{scenario} \noindentScenario \thescenario: $z_{i, \frac{j}{m},k}=\rho z_{i,\frac{j-1}{m},k}+\nu_{i,k}$, $k \in [K], j \in \{2,3,\ldots,d\}$, and $(\nu_{i,1}, \ldots, \nu_{i,K})^{\top} \sim_{\text {i.i.d}}$ student-$t$ with degree $2$ and its covariance matrix $\Sigma(1-\rho^2)$, with $\rho=0.9$ and $\rho_0 = 0.1$.

Furthermore, we examine how strong cross-covariate correlation affects the prediction and introduce the following scenario with an increased $\rho_0$.

\refstepcounter{scenario} Scenario \thescenario: $z_{i, \frac{j}{m},k}=\rho z_{i,\frac{j-1}{m},k}+\nu_{i,k}$, $k \in [K], j \in \{2,3,\ldots,d\}$, and $(\nu_{i,1}, \ldots, \nu_{i,K})^{\top} \sim_{\text {i.i.d }} \mathcal{N}\left(\boldsymbol{0}, \Sigma(1-\rho^2)\right)$ with $\rho=0.9$ and $\rho_0 = 0.9$. }

\textcolor{black}{It is important to note that these scenarios generate raw lagged covariates that enter directly into LASSO-UMIDAS, whereas LASSO-MIDAS and sg-LASSO-MIDAS require a MIDAS weighting transformation. Section (ref) discusses how this transformation affects these scenarios and provides a detailed comparison across them.}

To generate $T_i$, we first transform the covariates to their absolute values, which ensures that the distribution functions will be increasing in $t$ for all $\boldsymbol{Z}$. Then we let

for all $i\in[N]$, where $\mathbf{\zeta} \sim \text{Uniform}(0,1)$. The weighting schemes $\widetilde{\omega}_{k}(u), u \in [0,1]$ for $k =1, 2$ correspond to beta densities, respectively, equal to $\boldsymbol{\operatorname{Beta}}(1,3)$, $\boldsymbol{\operatorname{Beta}}(2,3)$, see ghysels2007midas,ghysels2019estimating,babii2022machine, for further details. This generation scheme guarantees that the survival function of $T$ satisfies (ref), where $ \boldsymbol{\theta}_0(t,s) $ is such that $ \boldsymbol{\theta}_{0,1}(t,s) =1 + \log(t-s) $, $\boldsymbol{\theta}_{0,1+j}(t,s) =(1+\log(t-s))\widetilde{\omega}_1\left(\frac{j-1}{d}\right),\ j\in[d]$, $\boldsymbol{\theta}_{0,1+d+j}(t,s) =(\log(t-s)-1)\widetilde{\omega}_2\left(\frac{j-1}{d}\right),\ j\in[d]$ and $\boldsymbol{\theta}_{0,k}(t,s) =0$ for all $k\in\{2d+2,\dots,K_z\}$. Remark that only the first two high-frequency covariates are relevant.

The censoring time $C_i$, for $i \in [N]$, is generated from a shifted exponential distribution. Specifically, we let $C_i = s + \text{Exp}_i$, where $\text{Exp}_i \sim \mathrm{Exp}(\gamma)$ is an exponential random variable with rate parameter $\gamma > 0$. The parameter $\gamma$ is chosen so that the resulting censoring rate, $(1/N)\sum_{i=1}^N \mathbbm{1}\{T_i > C_i\}$, is approximately $81\%$ in the simulated dataset, matching the censoring rate observed in the real dataset (see Section (ref)).

{\color{black} For the choice of the MIDAS weight function $W$ in the LASSO-MIDAS and sg-LASSO-MIDAS models, we use a dictionary of orthogonal polynomials shifted to the interval $[0,1]$, with size $L=3$, given by $W=\{w_1(u), w_2(u), w_3(u)\}$ for $u\in[0,1]$. The dictionary is constructed using a Gram–Schmidt orthogonalization of the power polynomials $\{1,x,x^2\}$ with respect to the measure $\mathrm{d}\mu(x)=(1-x)^{\alpha_{\text{poly}}}(1+x)^{\beta_{\text{poly}}}\mathrm{d}x$ on $[-1,1]$ where $\alpha_{\text{poly}} = \beta_{\text{poly}} =-1/2$. To obtain the basis on the unit interval, we apply the transformation $x=2u-1$, which yields the shifted polynomials $\{w_1(u), w_2(u), w_3(u)\}$. This construction is commonly referred to as a Gegenbauer polynomial dictionary babii2022machine in the MIDAS literature.\footnote{The function gb in the midasml R package computes an orthonormal polynomial basis for a given dictionary size and parameters $\alpha_{\text{poly}}$ and $\beta_{\text{poly}}$; see Online Appendix A of babii2022machine.} The use of orthogonal polynomial dictionaries helps reduce multicollinearity and improve numerical stability babii2022machine, though it may also change the variance of the MIDAS-weighted covariates, as discussed in Section (ref).

}

Regarding $t$, we set it to the following percentiles $t = \{t_1 = 10\%, t_2 = 30\%, t_3 = 50\%\}$ of the distribution of $\{T_i: T_i \text{ is uncensored}, i \in [N]\}$.

Concerning the evaluation of classification performance, Receiver Operating Characteristic (ROC) curves are widely used in the literature. However, traditional ROC curves are not fully suitable in this context due to censoring, where the status of firms is only partially observed. To address this limitation, we use the ROC curve estimator developed by heagerty2000time, which was specifically designed to evaluate classification performance effectively in the presence of censoring.

Evaluation metric: ROC curves with censoring

}

Recalling the definitions of sensitivity and specificity in the ROC curves, we see that in our model, both sensitivity, or the “true positive rate” (TPR), and specificity, or the “true negative rate” (TNR), are also functions that depend on $t$:

where $\Upsilon_i:= p\left(\widehat{\boldsymbol{\beta}}, \boldsymbol{X}_i\right) = \frac{\exp \left(\boldsymbol{X}_i^{\top}\widehat{\boldsymbol{\beta}}\right)}{1+\exp \left(\boldsymbol{X}_i^{\top}\widehat{\boldsymbol{\beta}}\right)}$ is the estimated probability.\footnote{When defining $\Upsilon_i$, we treat $\widehat{\beta}$ as fixed because it is estimated on the training set. The probabilities in $Se(c, t)$ and $Sp(c, t)$ are over the distribution of the test set.} The threshold $c$ is used to classify a firm as distressed if $\Upsilon_i > c $, or as non-distressed if $\Upsilon_i \leq c$, with $\mathbbm{1}\{T \leq t\}$ indicating whether the firm has failed by time $t$.

A ROC curve illustrates the full range of True Positive Rates (TPR) and False Positive Rates (FPR) across all possible threshold values $c$. A larger area under the ROC curve (AUC) signifies better performance in distinguishing between firms that have failed and those that have not. In practice, the status $\mathbbm{1}\{T \leq t\}$ in (ref) cannot be fully observed due to censoring. To address this issue, various ROC curve estimators have been proposed in heagerty2000time,cai2006sensitivity,heagerty2005survival,10.1093/biomet/asaa080.

Here, we employ the Nearest Neighbor estimator heagerty2000time to account for the censored data and evaluate the ROC curves. Let $$ \widehat{S}_{\kappa_N}(c, t)=\frac{1}{N} \sum_{i=1}^N \widehat{S}_{\kappa_N}\left(t \mid \Upsilon_i\right) \mathbbm{1}\{\Upsilon_i>c\}, $$ where $\widehat{S}_{\kappa_N}\left(t \mid \Upsilon_i\right)$ is a suitable estimator of the conditional survival function characterized by a parameter $\kappa_N$: $$ \widehat{S}_{\kappa_N}\left(t \mid \Upsilon_i\right)=\prod_{a \in \mathcal{T}_N, a \leq t}\left\{1-\frac{\sum_j \Psi_{\kappa_N}\left(\Upsilon_j, \Upsilon_i\right) \mathbbm{1}\{\widetilde{T}_j=a\} \delta_j}{\sum_j \Psi_{\kappa_N}\left(\Upsilon_j, \Upsilon_i\right) \mathbbm{1}\{\widetilde{T}_j \geq a\}}\right\}, $$ where $\mathcal{T}_N$ is a set of the unique values of $\widetilde{T}_i$ for observed events, $\delta_i = \mathbbm{1}\{T_i \leq C_i\}$ and $\Psi_{\kappa_N}\left(\Upsilon_j, \Upsilon_i\right)$ is a kernel function that depends on a smoothing parameter $\kappa_N$. Following the approach in heagerty2000time, we used a $0 / 1$ nearest neighbor kernel akritas1994nearest, $\Psi_{\kappa_N}\left(\Upsilon_j, \Upsilon_i\right)= \mathbbm{1}\{-\kappa_N<\widehat{F}_{\Upsilon}\left(\Upsilon_i\right)-\widehat{F}_{\Upsilon}\left(\Upsilon_j\right)<\kappa_N\}$ where $\widehat{F}_{\Upsilon}(\cdot)$ is the empirical distribution function of ${\Upsilon}$ and $2 \kappa_N \in(0,1)$ represents the percentage of individuals that are included in each neighborhood (boundaries). The resulting sensitivity and specificity are defined by: $$

$$ where $\widehat{S}_{\kappa_N}(t)=\widehat{S}_{\kappa_N}(-\infty, t)$. Both sensitivity and specificity above are monotone and bounded in $[0,1]$.

heagerty2000time used bootstrap resampling to estimate the confidence intervals for this ROC curve estimator. Motivated by the results of akritas1994nearest and cai2011robust, hung2010optimal discussed the asymptotic properties of the estimator and concluded that bootstrap resampling techniques can be used to estimate the variances of the proposed ROC curve. In practice, heagerty2000time suggested that the value for $\kappa_N$ is chosen to be $O\left(N^{-\frac{1}{3}}\right)$. In the present paper, we use the default value of the $\kappa_N$ produced in the documentation of the R package 'SurvivalROC', which is consistent with the choice found in blanche2013estimating. For further details on other ROC curve estimators in the survival analysis, we refer to kamarudin2017time.

{\color{black}

Discussion on the MIDAS-weighted covariates

In the simulation scenarios, increasing the autocorrelation strength $\rho$ of the original lags does not necessarily amplify the correlation within each group of MIDAS-weighted covariates. For example, although Scenario (ref) features stronger correlation among the original lagged covariates, the correlation matrix of the corresponding MIDAS--weighted covariates is nearly identical to that in Scenario (ref) (see Figure (ref)). This is because the orthogonal weighting functions in the MIDAS dictionary $W$ redistribute the dependence structure along orthogonal directions, so higher autocorrelation in the raw lags does not translate into stronger correlation among the transformed covariates.

Instead, the transformation could inflate the variance of these MIDAS-weighted covariates. To illustrate this, consider the $k$-th covariate with its $d$ high-frequency lags \[\widetilde{\boldsymbol{Z}}_{i,k} = \left(z_{i, s, k}, z_{i, s - \frac{1}{m}, k}, \ldots, z_{i, s - \frac{d-1}{m}, k}\right)^{\top},\] where each component has variance $1$ in both Scenarios (ref) and (ref). After multiplying by the MIDAS weighting matrix $W = \left(w_1, \ldots, w_L\right) \in \mathbb{R}^{d \times L}$, with $w_l \in \mathbb{R}^{d}$, we obtain the $k$-th group of MIDAS-weighted covariates $$W^{\top}\widetilde{\boldsymbol{Z}}_{i,k} = \left(\widetilde{\boldsymbol{Z}}_{i,k}^{\top}w_1, \ldots, \widetilde{\boldsymbol{Z}}_{i,k}^{\top}w_L \right)^{\top},$$ which is used in estimation. Unlike the original lags, the components of this transformed vector generally do not have variance equal to $1$. For example, $$Var\left(\widetilde{\boldsymbol{Z}}_{i,k}^{\top}w_1\right) = w_1^{\top}Cov\left(\widetilde{\boldsymbol{Z}}_{i,k} \right)w_1,$$ which depends on $W$ and the covariance structure of $\left(z_{i, s, k}, z_{i, s - \frac{1}{m}, k}, \ldots, z_{i, s - \frac{d-1}{m}, k}\right)^{\top}$, and hence varies with $\rho$ in the simulation design. In Scenarios (ref) and (ref), the stronger dependence among the original lagged covariates leads to larger variance in the corresponding MIDAS–weighted covariates.\footnote{\textcolor{black}{This pattern depends jointly on the covariance structure of $\widetilde{\boldsymbol{Z}}_{i,k}$ and the choice of the weighting matrix $W$. For instance, under a Toeplitz covariance structure combined with Chebyshev-based weights, specific patterns arise. In empirical applications, the true covariance structure of covariates is always unknown, so whether MIDAS necessarily increases their variance remains an open question.}} Table (ref) reports the average variance ratio of the $L$ elements in the first group ($k=1$) of MIDAS-weighted covariates between Scenario (ref) and Scenario (ref) based on the simulation, with similar patterns observed for other groups. Thus, a strong correlation among the raw lagged covariates increases the variance of the MIDAS-weighted covariates. This increase in the variance strengthens the signal relative to noise, thereby enhancing the predictive performance of MIDAS-based methods.

Notably, the MIDAS transformation is applied within each group, meaning to each original covariate and its lags, and does not change the correlation structure across different groups of covariates. Hence, we can evaluate how strong cross-covariate correlation affects the prediction of MIDAS-based methods by comparing Scenario (ref) and Scenario (ref). }

Prediction simulation results

We compute results for the three different LASSO-type regression methods. In the structured approach, sg-LASSO-MIDAS, each covariate and its high-frequency lags share the same group. Therefore, we have $K+1$ groups (one group corresponding to the intercept). Table (ref) presents the number of parameters (including the intercept) to be estimated in each of the three methods. It is evident that the two methods using MIDAS weights help mitigate the high-dimensional problem when $s \times m$ exceeds $L$.

We start by comparing the prediction results for sample sizes $N \in \{800, 1200\}$ across three simulation scenarios, followed by examining the recovery of the MIDAS weight function. To assess the prediction performance, we randomly split the simulated dataset into a training dataset ($80\%$) and a test dataset ($20\%$), ensuring that both sets maintain the same proportion of the event indicator $\delta_i(t) \mathbbm{1}{\{\widetilde{T}_i \leq t\}}$. We then calculate the estimated AUC in the test dataset, with the average estimated AUC obtained from $100$ simulated datasets for each sample size. The tuning parameters in the sg-LASSO-MIDAS and LASSO-MIDAS models are selected using $5$-fold stratified cross-validation to maximize the AUC on the validation fold, and the same procedure is applied to the LASSO-UMIDAS model. Specifically, we perform a grid search over the regularization parameter $\alpha$ in the sparse-group LASSO penalty, with values in the set $\{0, 0.1, 0.3, 0.5, 0.7, 0.9, 1\}$ and, as standard, $\lambda$ is chosen in a grid which follows liang_sparsegl.

Table (ref) reports the estimated average AUCs in the test dataset. {\color{black} To understand how different simulation scenarios affect the performance of a given method, we highlight two key factors. First, different simulation scenarios correspond to distinct models, as their true coefficients differ due to the prediction horizons $t$ being selected based on the generated data. For this reason, we report oracle AUCs for each scenario, computed by evaluating the AUC at the true data-generating parameters. The corresponding results have been added to Table (ref). Second, differences in covariate correlation strength across scenarios may affect the accuracy with which the true parameters are estimated. In what follows, we separately discuss the impact of covariate correlation on MIDAS-based methods (such as LASSO-MIDAS and sg-LASSO-MIDAS), drawing on the analysis in Section (ref), and on LASSO-UMIDAS.

We first focus on LASSO-MIDAS and sg-LASSO-MIDAS. As explained in Section (ref), it is expected that LASSO-MIDAS and sg-LASSO-MIDAS perform better in Scenario (ref) than in Scenario (ref). The oracle AUCs in Scenario (ref) are even higher than those in Scenario (ref) across different prediction horizons, which further explains the superior performance of LASSO-MIDAS and sg-LASSO-MIDAS. Similar patterns arise when comparing Scenario (ref) and Scenario (ref), with the only difference being that the covariates follow heavy-tailed distributions. As expected, Scenario (ref) features stronger correlations across groups of MIDAS-weighted covariates, leading to noticeably weaker performance of LASSO-MIDAS and sg-LASSO-MIDAS compared with Scenario (ref). The lower oracle AUCs in Scenario (ref) further reinforce this pattern. Although Scenarios (ref) and (ref) exhibit the highest oracle AUCs overall, LASSO-MIDAS and sg-LASSO-MIDAS perform worse there than in Scenarios (ref) and (ref) due to the heavy-tailed covariates.

As for LASSO-UMIDAS, it exhibits performance patterns similar to those of LASSO-MIDAS and sg-LASSO-MIDAS across the scenarios. It attains a higher AUC in Scenario (ref) than in Scenario (ref), partly reflecting the higher oracle AUCs. This is intuitive: Scenario (ref) features stronger within-group correlations among the active covariates but similarly low cross-group correlations, which allows LASSO-UMIDAS to approximate the predictive signal by selecting a few representative covariates in the active groups and therefore yields more stable predictions. In contrast, under the weaker within-group correlation of Scenario (ref), the signal is more diffusely distributed across the active covariates, leading LASSO-UMIDAS to select too few predictors and produce less stable predictions with lower AUC. To quantify differences in signal strength across scenarios, Table (ref) reports the relative signal-to-noise ratios (SNRs) for LASSO-UMIDAS (shown for \(N = 800\), with similar results for \(N = 1200\)). We measure the SNR as $3\operatorname{Var}\!\left(\boldsymbol{Z}^{\top}\boldsymbol{\theta}_0\right)/\pi^2 $ using the standard normalization that the logistic error has variance $\pi^2/3$. As noted by chardon2024finite, stronger signal strength facilitates classification. The substantially higher SNRs in Scenario (ref) explain why LASSO-UMIDAS achieves more accurate predictions there compared with Scenario (ref). Although we do not report SNRs for Scenarios (ref) and (ref), the same qualitative pattern is observed: these scenarios differ from Scenarios (ref) and (ref) only in using heavy-tailed covariates, and the resulting predictions follow a similar trend. Table (ref) further supports that LASSO-UMIDAS performs better in Scenario (ref) than in Scenario (ref), as the SNRs are higher in Scenario (ref), where the strong correlation among active lags appears to be particularly beneficial. Although high within-group correlation among active covariates tends to improve prediction, strong cross-group correlation between active and inactive groups can diminish predictive performance. The decrease in AUCs from Scenario (ref) to Scenario (ref) reflects the effect of stronger cross-group dependence. When comparing Scenarios (ref) and (ref) with Scenarios (ref) and (ref), the degraded performance in the latter can be attributed to the heavy-tailed covariates. For $N = 800$, LASSO-UMIDAS performs similarly in Scenarios (ref) and (ref), but the AUCs across all three prediction horizons are consistently slightly closer to the oracle in Scenario (ref), indicating a modest advantage in this scenario. LASSO-UMIDAS also consistently outperforms in Scenario (ref) compared with Scenario (ref), as expected. When the sample size increases to $N = 1200$, the same pattern persists: LASSO-UMIDAS shows better performance in Scenarios (ref) and (ref), while Scenarios (ref) and (ref) remain affected by heavy-tailed covariates.

For general context, as discussed in hebiri2012correlations,lassoprediction, the predictive performance of Lasso-type methods can remain stable under various levels of correlation provided that the penalty parameter is appropriately tuned, for example, through cross-validation. Moreover, since the five scenarios correspond to different underlying data-generating models, our primary objective is to compare how the MIDAS-based methods perform relative to LASSO-UMIDAS within each scenario.

We therefore focus on method comparisons within the same scenario. As shown, sg-LASSO-MIDAS achieves the highest AUCs across different simulation scenarios. Both sg-LASSO-MIDAS and LASSO-MIDAS, using weight function approximations, outperform LASSO-UMIDAS. LASSO-UMIDAS generally demonstrates the poorest predictive performance across all scenarios. As expected, the predictive performance improves with an increase in sample size $N$. These results remain robust as the autocorrelation of original lagged covariates $\rho$ increases from $0.1$ to $0.9$ or the cross-covariate correlation $\rho_0$ increases from $0.1$ to $0.9$. Although all three methods perform less effectively with heavy-tailed covariates, sg-LASSO-MIDAS continues to outperform the others. In Tables (ref), (ref), (ref), (ref), and (ref) of Online Appendix (ref), we report additional results for the estimation accuracy of the true parameters. } It is worth noting that the increase of the parameter estimation accuracy with the sample size is not particularly large, as the high censoring rate in the simulated datasets limits the average increase in the number of uncensored firms ($\mathbbm{1}\{T_i \leq C_i\} = 1$) to only about $76$ as the sample size $N$ grows from $800$ to $1200$.\footnote{In simulation results not shown for brevity, when the dataset's censoring rate is approximately $30\%$, notable improvements both in parameter estimation and estimated AUC are observed as the sample size $N$ increases from $800$ to $1200$.}

{\color{black} In addition to the prediction performance, we examine the variable selection performance of sg-LASSO-MIDAS across different scenarios.

Recall that we have $K = 50$ covariate groups, with the first two groups being active and the remaining $48$ groups being inactive. We define a group as selected if any estimated coefficient in the group is nonzero. That is, for group $k$: \[ A_k =

. \] The performance metrics for variable selection are defined as follows:

Table (ref) reports the TPR for sg-LASSO-MIDAS across different scenarios. The pattern of variable selection is similar to the trends observed in the estimated AUCs in Table (ref). As expected, sg-LASSO-MIDAS achieves its best TPR in Scenario (ref), where high within-group correlation is effectively handled, and the SNR is higher as a result of the MIDAS transformation relative to the other scenarios. In contrast, variable selection deteriorates in Scenarios (ref), (ref), and (ref), due to the presence of heavy-tailed covariates and high cross-covariate correlations. Across all scenarios, performance improves with larger sample sizes, as increased data stabilizes estimation and variable selection.

Overall, the simulation evidence strongly supports the advantage of using MIDAS weighting and incorporating the internal structure of covariates in high-dimensional settings. }

{\color{black}

Inference simulation results

To assess the inference performance, we keep the data generation process in Section (ref), set the first group of MIDAS-weighted parameters as of interest, and test

We implement a Wald test, using a nominal significance level of $5\%$ in this simulation. We scale the Beta density function $\widetilde{\omega}_{1}(u)$ by multiplying it with a constant $a \in \{ 0, 0.1 \}$. The tuning parameters for both sg-LASSO-MIDAS and the nodewise LASSO are selected using $5$-fold cross-validation to maximize the likelihood score on the validation folds. For the regularization parameter $\alpha$ in the sparse-group LASSO, we use a relatively small grid set $\{0, 0.5, 1\}$ to alleviate computational burden. The sample size is $N = 1200$ and the simulation is run $500$ times. Unlike for the prediction exercise, we do not split the data here and use the full dataset to estimate. All other settings follow those in Section (ref).

As shown in Table (ref), the proposed de-sparsified sg-LASSO-MIDAS attains an empirical test size close to the nominal $5\%$ level in both Scenario (ref) and Scenario (ref). The empirical power in Scenario (ref) is higher than in Scenario (ref), which is consistent with the larger signal-to-noise ratio discussed in Section (ref). In Scenarios (ref) and (ref), which involve heavy-tailed covariates, a similar pattern holds: the method performs better in terms of both size and power in Scenario (ref), which has a larger signal-to-noise ratio. Empirical power drops markedly when cross-covariate correlation increases, as in the comparison of Scenario (ref) and Scenario (ref). Comparing Scenario (ref) with Scenario (ref), the latter shows higher empirical power but also a substantial deviation of empirical size from the nominal $5\%$ level, indicating notable size inflation. Both Scenario (ref) and Scenario (ref) maintain empirical sizes close to the nominal level, with Scenario (ref) achieving higher empirical power.

}

Empirical application

Data

We construct a dataset of all publicly traded Chinese manufacturing firms listed on the Shanghai and Shenzhen Stock Exchanges. These firms' financial statuses are classified as either Special Treatment (ST) or No-ST.\footnote{The initial public offering (IPO) dates of these firms fall between $1985$, January $1^{\text{st}}$ and $2015$, December $31^{\text{st}}$.} A firm is designated as an ST firm if it meets any of the following criteria: i) two consecutive years of earnings are negative; ii) one recent year of earnings is negative and the most recent year of equity is negative; iii) the most recent year’s audited financial statements conclude with substantial doubt; and iv) other situations identified by the stock exchange as abnormal activities or a high risk of delisting. According to li2021chinese, ST status is a reliable indicator of financial distress in China. Therefore, we use the ST indicator as a proxy for a firm's financial distress.

The dataset is sourced from the IFIND database \url{https://www.hithink.com/ifind.html}, one of China's leading financial data providers. The database contains mostly manually extracted data, covering financial data such as stocks, bonds, funds, futures, and indexes. \textcolor{black}{Detailed information about the dataset can be found in Online Appendix (ref).} Additionally, we have developed an R package, Survivalml, which is publicly available at \url{https://github.com/Wei-M-Wei/Survivalml}, and its instructions are provided as well.

The raw dataset consists of $1614$ companies, of which $299$ were classified as ST and $1315$ as No-ST. The data cover the period from January $1^{\text{st}}, 1985$ to December $31^{\text{st}}, 2020$, with firms entering the sample from their respective IPO dates within this window.\footnote{There are no mergers in the dataset.} The dataset exhibits a censoring rate of approximately $81\%$. We collect $57$ quarterly measured financial variables, categorized into $8$ types (number of covariates in each type), as follows: Operation-Related $(6)$, Debt-Related $(10)$, Profit-Related $(16)$, Potential-Related $(6)$, Z-score Related $(5)$ altman1968financial, Capital-Related $(6)$, Stock-Related $(5)$, and Cash-Related $(3)$. Table (ref) provides detailed information on these financial variables; see Online Appendix (ref) for further details. Figure (ref) presents the distribution of IPO, first-time-to-be ST, and censored firms across different years of the raw dataset. Many of them were listed in $2010$ and $2011$, several firms were publicly listed in $2013$, and the financial distress firms seem to be distributed evenly between $1999$ and $2020$.\footnote{Since China put froze IPOs in $2013$, there were only a limited number of IPO firms in this year.}

We construct the sub-dataset in which all firms have survived for at least $s$ years. The goal is to use these $s$ years of information to predict whether a firm will fail within $t$ years.

Estimation procedure

We now describe the estimation procedure in the empirical application. First, we note that all public firms report their financial information with a one-quarter delay. Consequently, if a firm has survived for $s$ years, only $s \times 4 - 1$ quarters' worth of financial information will be available for analysis.

Let $z_{i, s - \frac{j-1}{m}, k}$ represent the $k$-th financial covariate of firm $i$, measured at time $s - \frac{j-1}{m}$, where $j = 2, \ldots, d$, and $d = s \times m$. We organize all the lags of the covariate into a group vector $\widetilde{\boldsymbol{Z}}_{i,k}$: $$ \widetilde{\boldsymbol{Z}}_{i,k} = \left(z_{i, s - \frac{1}{m}, k}, z_{i, s - \frac{2}{m}, k}, \ldots, z_{i, \frac{1}{m} ,k}\right)^{\top}, \quad i \in [N], \quad k \in [K], $$ where $z_{i,\frac{1}{m}, k}$ refers to the $k$-th covariate measured in the next quarter following the firm's IPO date, and $m = 4$ denotes the quarterly frequency of the financial covariates.

Next, we aggregate the lagged covariate vector $\widetilde{\boldsymbol{Z}}_{i,k}$ using a dictionary $W$, which consists of Gegenbauer polynomials shifted to the interval $[0,1]$ with parameter $\alpha_{\text{poly}} = \beta_{\text{poly}} = -\frac{1}{2}$ and size $L = 3$.\footnote{These polynomials are also known as Chebyshev polynomials on the interval $[0,1]$.} This choice of $W$ coincides with the specification used in Section (ref).

Finally, we construct the covariate matrix $\boldsymbol{X}$ as follows: $$\boldsymbol{X} =\left(\boldsymbol{X}_{1} , \boldsymbol{X}_{2}, \ldots, \boldsymbol{X}_{N}\right)^{\top},$$ where each $\boldsymbol{X}_i$ $=\left(1, \widetilde{\boldsymbol{Z}}_{i,1}^{\top} W, \widetilde{\boldsymbol{Z}}_{i,2}^{\top} W, \ldots, \widetilde{\boldsymbol{Z}}_{i,K}^{\top} W\right)^{\top}, i \in [N]$. This matrix $\boldsymbol{X}$ is then used in sg-LASSO-MIDAS and LASSO-MIDAS. Notice that we include the intercept term but do not penalize it in the estimation procedure.

We compare the performance of firm distress predictions using the following methods.

\noindentLogistic regression. \textcolor{black}{As a benchmark, we consider a simple unpenalized logistic regression using only the lags of the first Z-score covariate. We solve the empirical version of (ref), in which the function $H$ is estimated using the Kaplan-Meier estimator. This is considered a reasonable starting point for distress predictions. The total number of parameters, including the intercept, to estimate is $1 + (d - 1)$, where $d-1$ is the number of available lags for each covariate.}

\noindentLASSO-U (LASSO-UMIDAS). We estimate $ d - 1 = s \times m -1$ coefficients per group covariate $\widetilde{\boldsymbol{Z}}_{i,k} \in \mathbb{R}^{d-1}, k \in [K]$, using the unstructured LASSO estimator. The total number of parameters to estimate is $1 + K \times (d - 1)$, where $s$ is the number of years survived by the firm, and $m$ represents the annual sampling frequency of each covariate.

\noindentLASSO-M (LASSO-MIDAS). Each high-frequency covariate and its $d-1$ lags are grouped into $\widetilde{\boldsymbol{Z}}_{i,k} \in \mathbb{R}^{d-1}, k \in [K]$. We aggregate the group covariate $\widetilde{\boldsymbol{Z}}_{i,k}$ using Gegenbauer polynomials $W \in \mathbb{R}^{(d-1) \times L}$. We apply a Lasso penalty to induce sparsity. The total number of parameters to estimate is $1 + K \times L$, where $L$ is the size of the Gegenbauer polynomial dictionary.

\noindentsg-LASSO-M (sg-LASSO-MIDAS). Similarly to LASSO-MIDAS, each high-frequency covariate and its $d-1$ lags form a group $\widetilde{\boldsymbol{Z}}_{i,k} \in \mathbb{R}^{d-1}, k \in [K]$, which is aggregated using Gegenbauer polynomials $W \in \mathbb{R}^{(d-1) \times L}$. Instead of using a Lasso penalty, we use the sparse-group Lasso penalty to induce sparsity in the group covariates. \textcolor{black}{The group structure matches that used in the simulation, where each group corresponds to the MIDAS-weighted covariates. Alternative group structures are also examined, and the corresponding results are reported in Online Appendix (ref).} The total number of parameters to estimate is $1 + K \times L$.

The choice of $s$ dictates the historical information captured in the covariate matrix $\boldsymbol{X}$, while $t$ denotes the prediction horizon. The existing literature on firm distress prediction, particularly in the United States, often examines prediction period $t-s$ ranging from $1$ quarter to $2$ years cole2012deja. In practical applications, such as for bank regulators, models need to identify potential failures well in advance. For example, audrino2019predicting developed a MIDAS-type method with prediction periods $t-s$ equal to $1$ and $2$ years.

For our empirical application, we select a reference period of $s = 6$ years. Using the firm classification criteria for Special Treatment outlined in Section (ref), we establish prediction horizons of $t = 8, 8.5, 9$ years to forecast firm distress within these intervals. To investigate longer forecast periods, we also consider an additional case with $s = 10$ years and prediction horizons of $t = 13, 13.5, 14$ years. Longer prediction horizons provide information on the risks of long-term financial distress. As highlighted by li2021chinese, these prediction horizons are critical for accurately forecasting firm financial distress in China. They also offer meaningful and practical benchmarks for evaluating firm failure prediction models.

In practice, missing data in financial variables can arise due to various factors, including inconsistent reporting practices between firms, differing regulatory requirements, incomplete disclosures, and delays in data availability after IPOs. Given the substantial amount of missing data in the raw dataset, we construct a complete sub-dataset for each $s$ by selecting firms with consistent $s$-year observations. While common approaches for handling missing data, such as removing variables or firms with missing values, are widely used, these methods often result in retaining too few firms or variables for meaningful analysis. Furthermore, the firms with observable status in the sub-dataset play a critical role in the predictive modeling process. To address these challenges, we propose an algorithm that balances dimensionality and the number of uncensored firms in the selected sub-dataset, as outlined in Algorithm (ref) of Online Appendix (ref).

To evaluate the prediction performance of the different methods, we randomly split the dataset into in-sample ($80\%$) and out-of-sample ($20\%$) sets, ensuring that both sets maintain the same proportion of the event indicator $\delta_i(t) \mathbbm{1}\{\widetilde{T}_i \leq t\}$. The tuning parameters for sg-LASSO-MIDAS, LASSO-MIDAS, and LASSO-UMIDAS are selected using stratified $5$-fold cross-validation, where the optimal parameters are those that maximize the AUC on the validation fold.\footnote{In this paper, unless specified otherwise, the default choice for cross-validation is to maximize the AUC on the validation fold.} Additionally, as an alternative, cross-validation to maximize the likelihood score is also investigated. The AUC estimator employed in this procedure follows the method described in Section (ref)}. Specifically, we perform a grid search over the regularization parameter $\alpha$ in the sparse-group LASSO penalty, with values in the set $\{0, 0.1, 0.3, 0.5, 0.7, 0.9, 1\}$ and, as standard, $\lambda$ is chosen in a grid which follows liang_sparsegl. This process is repeated $10$ times as a robustness check, each time using a different random split of the data. All models are trained on the same training set and evaluated on the same test set.

For each split, the AUC is computed on the out-of-sample data, and the out-of-sample data is then bootstrapped $1000$ times to calculate the AUC for each bootstrap sample. The AUC values for each bootstrap sample are subsequently averaged across the $10$ different splits, resulting in $1000$ averaged AUC values. The final performance is reported as the overall average AUC, along with a two-side $95\%$ confidence interval, which is calculated based on these $1000$ bootstrapped averages. This approach ensures a robust performance evaluation by accounting for variability in the data and model performance.\footnote{We note that the bootstrap approach is not theoretically validated for the regularized estimators we consider in this paper.}

On top of the simple logistic regression and the LASSO-UMIDAS, LASSO-MIDAS, and sg-LASSO-MIDAS with cross-validation for the AUC or the likelihood score, we consider other alternative approaches.

\paragraph{Macro data augmented prediction.} We first assess whether incorporating macroeconomic data can enhance the accuracy of distress prediction models. The macroeconomic dataset for China is sourced from the Federal Reserve Bank of Atlanta’s China Macroeconomy Project (\url{https://www.atlantafed.org/cqer/research/china-macroeconomy#Tab2}), which provides a comprehensive set of macroeconomic variables relevant to the Chinese economy. The dataset includes $98$ macroeconomic variables, measured quarterly, and spans the same time period as the financial data collected for the firms in our study.

To merge the macroeconomic data with the financial dataset, we select only those macroeconomic variables that do not have missing values across all firms within each financial sub-dataset. Since the sub-datasets differ based on the value of $s$, the set of macroeconomic variables selected will vary accordingly for each sub-dataset. Furthermore, we use the same MIDAS dictionary $W$ for the macroeconomic covariates as for the financial covariates, ensuring consistency in the aggregation of high-frequency data over time.

Table (ref) summarizes the details of the two sub-datasets categorized by different values of $s$. For the sub-dataset with $s=6$ years, we use all available information across each firm’s entire survival period, allowing us to leverage the maximum historical data available for firms with $6$ years of survival. In contrast, for the sub-dataset with $s=10$ years, we restrict the covariates to those from the last $4$ years of each firm’s survival period. This adjustment is necessary because firms that have survived for $10$ years were generally listed in the $1990$s, and significant missing data is often observed in the early years following their IPOs. By focusing on the most recent $4$ years, we ensure better data quality and a more robust analysis.

\paragraph{Oversampling.} In addition, we apply an oversampling technique to address the imbalance in the dataset caused by the high censoring rate, which results in an unequal proportion of firms experiencing distress versus those that are not. This imbalance could adversely affect the performance of distress prediction models, as the minority class (distressed firms) may be underrepresented. Since the empirical dataset has a high censoring rate, we face a class imbalance between those firms that eventually experience distress $ \mathbbm{1}\{T_i \leq C_i\}\mathbbm{1}\{\widetilde{T}_i \leq t\} = 1 $ and those that do not or we do not observe $\mathbbm{1}\{T_i \leq C_i\}\mathbbm{1}\{\widetilde{T}_i \leq t\} = 0$. To balance this, for the training dataset, we randomly duplicate the observations from the minority class (firms that experience distress) until the proportion of distressed firms reaches $15\%$ of the training dataset. This step helps mitigate the imbalance and ensures that the model is exposed to a sufficient number of distressed firms during training. Tuning parameters are selected using $5$-fold stratified cross-validation, where the optimal parameters maximize the likelihood score.

\paragraph{Does censoring matter for prediction?} We compare with an approach that applies LASSO-UMIDAS, LASSO-MIDAS, and sg-LASSO-MIDAS to the sub-dataset where censored firms with censoring time smaller than $t$ ($\mathbbm{1}\{ C_i < T_i\}\mathbbm{1}\{ C_i < t\}=1$) have been removed.\footnote{Given the prediction horizon $t$, the distress status of censored firms cannot be observed if their censoring time is shorter than $t$ and, clearly, censored firms with $\mathbbm{1}\{C_i \geq t\}$ are not distressed.} This is the approach usually taken in the literature, since it allows for ignoring censoring, see the discussion in the literature review of the introduction. The limitation of this procedure is that it does not use all observations, resulting in a loss of precision.

Application results

The results are presented in Tables (ref) and (ref). The LASSO-MIDAS and sg-LASSO-MIDAS consistently outperform the LASSO-UMIDAS, which aligns with our expectations and the logistic regression benchmark. For $s=6$ years, sg-LASSO-MIDAS and LASSO-MIDAS perform essentially equally, with only minor differences across prediction horizons. When $s=10$ years, however, sg-LASSO-MIDAS exhibits a slight but consistent performance advantage over LASSO-MIDAS, indicating that the sparse-group Lasso regularization is particularly beneficial when incorporating a larger historical window of data.

When we compare the performance of models based on cross-validation using different metrics, we observe that cross-validation based on the AUC generally yields better results than cross-validation based on likelihood scores. This is logical since our target measure is the AUC itself.

Additionally, while integrating macroeconomic data does not improve prediction performance over the purely financial model when $s = 6$ years, it enhances performance when $s = 10$ years. This suggests that macroeconomic variables become more relevant with a larger historical window, offering supplementary information that helps improve prediction accuracy, especially for firms with longer survival periods. However, oversampling does not seem to provide any additional benefit in improving prediction performance.

When we remove censored firms with $C_i < t$, the performance of our methodologies deteriorates across all scenarios, emphasizing the importance of properly accounting for censoring in predictive modeling.

\textcolor{black}{We also evaluate how different covariate group structures influence the performance of sg-LASSO-MIDAS, with results reported in Tables (ref) and (ref) of Online Appendix (ref). The performance of sg-LASSO-MIDAS is generally robust to variations in group structure, consistently outperforming LASSO-UMIDAS and, under certain group specifications, surpassing LASSO-MIDAS more frequently than under the baseline grouping; see Online Appendix (ref) for a description of the group structures and further discussion.}

To further assess the performance difference, we conduct a pairwise comparison test, a widely used method for comparing two AUCs. Slightly modifying the approach in han, robin_proc_2011, we specifically test whether the estimated AUC of sg-LASSO-MIDAS is no larger than \textcolor{black}{those of LASSO-MIDAS and LASSO-UMIDAS.}\footnote{We have $1000$ bootstrapped average AUCs for each method as described before. The p-value is calculated as the proportion of sg-LASSO-MIDAS's AUC values that are smaller than those of another method.} Table (ref) shows that the improvement of the sg-LASSO-MIDAS over the LASSO-UMIDAS is statistically significant at least at the $10\%$ significance level across all scenarios, with the largest gap observed when $s = 10$ and $t = 13$ years. \textcolor{black}{As for the comparison between sg-LASSO-MIDAS and LASSO-MIDAS, no significant difference is found when $s = 6$ years. However, when $s =10$ years, sg-LASSO-MIDAS significantly outperforms LASSO-MIDAS at the $10\%$ significance level for prediction horizons $ t = 13,13.5$ years.} We also conduct a pairwise comparison between the sg-LASSO-MIDAS applied to the dataset with and without censored firms satisfying $C_i < t$. For the scenarios where $s = 6$, $t = 8.5$ years, and $s = 10$, $t = 13$ years, sg-LASSO-MIDAS performs better on the full dataset than on the dataset without censored firms satisfying $C_i < t$, though the difference is not statistically significant. However, in other scenarios, including censoring significantly improves model performance, with results being statistically significant, at least at the $5\%$ significance level. These findings strongly support the advantages of using MIDAS weights, considering the group structure of covariates, and incorporating the censoring information in practice. Overall, the empirical results highlight the superiority of the sg-LASSO-MIDAS across different scenarios.

To better understand which covariates are useful for prediction, we examine the financial types selected by the sg-LASSO-MIDAS, as illustrated in Figure (ref) of Online Appendix (ref). Financial variables related to the $Z$-score appear to play a pivotal role across all prediction horizons in forecasting firm distress. This observation aligns with prior research altman1968financial, as the $Z$-score model has been widely employed in both academic studies and industry to predict corporate defaults altman2017financial. Further details on the selected financial covariates are presented in Figures (ref) and (ref) of Online Appendix (ref).

{\color{black} Furthermore, following Section (ref), we separately test each financial covariate for statistical significance at the $5\%$ level in distress prediction. Specifically, for each covariate, we conduct a Wald test on its associated vector of MIDAS coefficients, estimated using the de-sparsified sg-LASSO-MIDAS on the full empirical dataset, as in Section (ref). All the MIDAS settings remain the same as those in Section (ref). For the regularization parameter $\alpha$ in the sparse-group LASSO, we consider the set $\{0, 0.1, 0.3, 0.5, 0.7, 0.9, 1\}$. Tuning parameters for sg-LASSO-MIDAS and the nodewise LASSO are chosen via $5$-fold stratified cross-validation to maximize the validation likelihood. As shown in Figures (ref) and (ref) of Online Appendix (ref), statistically significant covariates are marked in blue for each prediction horizon. Consistent with the literature altman1968financial, altman2017financial, we find that the five Z-score covariates (denoted by $X_1,\ldots,X_5$) play prominent roles, while several other financial covariates also have a notable influence. When $s = 6$ years, $X_2$ (Retained Earnings\,/\,Total Assets) and $X_4$ (Market Value of Equity\,/\,Total Liabilities) are statistically significant across all prediction horizons, whereas $X_1$ (Working Capital\,/\,Total Assets) is significant at only one prediction horizon. Beyond the Z-score covariates, several other covariates are also significant. For example, Earnings Before Interest, Taxes, Depreciation, and Amortization\,/\,Total Liabilities is consistently significant across all horizons. For $s = 10$ years, $X_4$ remains significant across all horizons, while $X_2$ and $X_3$ (Earnings Before Interest and Taxes\,/\,Total Assets) are significant at only one prediction horizon. Other financial covariates, such as Net Profit on Assets and the \emph{Interest-Bearing Debt Ratio}, are also influential at most prediction horizons. }

Conclusion

This paper presents a novel approach to corporate survival analysis, addressing the challenges of high-dimensional censored data sampled at both consistent and mixed frequencies.

The first major contribution is the introduction of the sparse-group LASSO estimator for high-dimensional censored MIDAS logistic regressions, \textcolor{black}{along with a corresponding de-sparsified estimator for inference.} The proposed framework effectively accommodates hierarchical data structures and facilitates variable selection both within and across groups, unifying classical LASSO and group LASSO.

Secondly, we develop the theory for logistic regression with high-dimensional censored data sampled at different frequencies. To extend the existing literature with assumptions on fixed design or isotropic conditions of the sub-Gaussian covariates, we develop the non-asymptotic properties of the proposed sparse-group LASSO estimator for censored, heavy-tailed data. The framework is readily extendable to generalized linear models with structured sparsity estimators. We also explicitly consider the approximation error, which, to the best of our knowledge, is a novel contribution in the context of logistic regression. This error may arise from various sources, including approximations in the MIDAS weight function and/or deviations from exact sparsity. \textcolor{black}{In addition to the penalized estimator, we analyze the corresponding de-sparsified sparse-group LASSO estimator, showing that it is asymptotically unbiased while properly accounting for censoring, MIDAS approximation error, and heavy-tailed covariates, and that its asymptotic variance is affected by censoring, which is not previously studied in the literature.}

A key practical contribution is an application to a comprehensive dataset of publicly traded Chinese manufacturing firms, integrating survival and censoring time information alongside numerous high-frequency financial covariates. Empirical findings show that sg-LASSO-MIDAS consistently outperforms unstructured LASSO methods across various scenarios. Notably, the inclusion of censoring information significantly enhances prediction performance, providing valuable insights for predicting firm distress under real-world conditions. \textcolor{black}{Beyond prediction, the de-sparsified sg-LASSO-MIDAS further enables valid statistical inference, allowing us to identify financial covariates that are statistically significant predictors of firm distress and thereby providing practical guidance for decision makers.}

Overall, the methodologies developed in this paper have broad applicability beyond corporate distress prediction. The integration of logistic models, MIDAS, and regularized machine learning techniques holds promise for applications in areas such as disease diagnosis, solvency evaluation, fraud detection, customer churn analysis, and labor market studies.

Acknowledgments

The authors thank the Associate Editor and two anonymous referees, Christophe Croux, Geert Dhaene, Daniel Gutknecht, Onno Kleen and Yoshimasa Uematsu for helpful comments as well as seminar participants at Tilburg University and Copenhagen Business School, and conference attendants at the International Association for Applied Econometrics 2025, the Leuven Statistics Days 2025, the International Conference on Econometrics and Statistics 2025, and Financial Econometrics Meets Machine Learning 2025. The authors are listed in order of contribution.

{

Funding sources

Wei Miao gratefully acknowledges financial support from the China Scholarship Council through the grant 202306130036. Jad Beyhum gratefully acknowledges financial support from the Research Fund KU Leuven through the grant STG/23/014. Jonas Striaukas gratefully acknowledges the financial support from the European Commission, MSCA-2022-PF Individual Fellowship, Project 101103508. Project Acronym: MACROML. Ingrid Van Keilegom acknowledges support from the FWO (Research Foundation Flanders) through the projects G0I3422N and G047524N. }

Supplementary material

\paragraph{Online Appendix:} \textcolor{black}{Proof of Theorem (ref), Lemma (ref), and Theorem (ref), empirical dataset pre-processing algorithm, additional empirical and simulation results, and details on the empirical dataset (.pdf file).} \paragraph{R package and Replication files:} \textcolor{black}{R package ‘Survivalml’ that implements our method, and all simulation and application code are available on \url{https://github.com/Wei-M-Wei/Survivalml}.}