Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
70,537 characters · 15 sections · 99 citation commands
Free Discontinuity Regression With an Application to the Economic Effects of Internet Shutdowns
\singlespacing
\thispagestyle{empty}
\thispagestyle{empty} \setcounter{page}{1}
\onehalfspacing
In many settings of interest, discontinuous changes in regression surfaces can reveal fundamental insights into the underlying problem structure, with both location and size of discontinuities typically unknown and dependent on multiple variables. Such discontinuities arise in many contexts, including algorithmic decision-making in finance and education, marketing segmentation, and tipping points in climate, urban sorting, financial markets, and spatial gene-expression profiling argyle2020monthly, brunner2021effects, kuo2002integration, li2010contextual, lamberson2012tipping, card2008tipping, scheffer2009early, hansen2017regression, caetano2017school, aue2009break, chen2023scs, hou2023statistical.
We introduce Free Discontinuity Regression (FDR), a method that reconstructs the regression surface and the exact geometry and magnitude of its jumps in one shot. FDR segments the domain into smooth and discontinuous regions without any point-wise smoothness assumption. In contrast to wavelet smoothing, fused-lasso variants, or gradient-threshold rules, our approach embeds smoothing and boundary detection in the same optimization, so the jump set and its sizes emerge as primary estimands rather than as results of a second-stage thresholding procedure, allowing us to develop the first formal identification guarantees for both jump location and size estimation in any dimension.
Our starting point is the Mumford-Shah (MS) functional mumford1989optimal, a workhorse model from mathematical computer vision and extension of classic discrete approaches geman1984stochastic, which jointly segments and denoises while incorporating spatial regularity. The MS functional arises naturally as a solution to the limitations of simpler edge detection and thresholding methods canny1986computational, kass1988snakes, which almost surely misclassify discontinuities in noisy settings chan2005image and typically require multiple estimation steps, complicating statistical inference.\footnote{For empirical validation of the MS functional's superiority over simpler approaches, see wang2012robust, chan2001active, lucas2022hyperparameter, strekalovskiy2012convex.} The MS functional is non-convex, so its critical points need not be global optima. This poses a threat to the reproducibility of estimates and the ability to characterize convergence to the true population function. To resolve it, we employ the method of calibrations alberti2003calibration to obtain a convex relaxation with global solutions, implementable via primal-dual algorithms chambolle2011first.
Most classical approaches to the MS functional treat the noise as already baked into fixed image $f$. Although $f$ may deviate from the true underlying signal $u$, the analysis typically does not place a probability distribution over all possible noisy images; rather, $f$ is taken as a single, deterministic function. In contrast, our statistical viewpoint treats the observed function $Y_i$ as doubly random: the function values are observed with random noise, and live in a random point cloud. caroccia2020mumford recently considered a similar setting for the Mumford-Shah functional on random graphs generated by random point clouds, where the noise is additive and i.i.d. In contrast, we work with a random grid on a random point cloud, and allow for the noise to be correlated. This is important in spatial settings like our application, and even in the classical image setting, where we can expect noise in blurry segments of the image to be correlated. Moreover, we provide sharp identification results in any dimension, in addition to $\Gamma$-convergence---the first results of their kind, to our knowledge. These complement and extend recent results for discrete approximations of the classical non-convex MS functional caroccia2020mumford, ruf2019discrete, morini2002global, richardson1992limit and related functionals chambolle2021learning, garcia2016continuum, trillos2017new, hutter2016optimal, hu2022voronoigram.
Figure (ref) demonstrates our method on 1D-3D simulated examples, showing accurate recovery of both functions and discontinuities from noisy data. We note that, while our main contribution lies in estimating multivariate discontinuities, our method also performs remarkably well in the one-dimensional case, handling multiple jumps of different sizes effortlessly without requiring user input on the number of jumps or the minimum jump size. Further below, we apply the method to estimate the effects of an internet shutdown in Rajasthan, India, in 2021. This one-day shutdown, implemented to prevent exam cheating, created sharp geographic discontinuities in connectivity. Using anonymized mobile activity around economic points of interest and cached location data from satellite-connected apps, we estimate a 25-35% reduction in economic activity along the shutdown boundary. These findings substantially exceed previous estimates and highlight the asymmetry between internet expansion and shutdowns in modern economies, as well as the critical importance of reliable digital infrastructure for economic performance.
Our method advances several literatures. For a concise overview of how it compares to the main competing methods, see Table (ref). The rich literature on change point detection methods typically focuses on discontinuity locations rather than magnitudes, or on the one-dimensional setting page1954continuous, muller1992change, braun2000multiple, killick2012optimal, porter2015regression, donoho1994ideal, harchaoui2008kernel, while structural break estimation primarily addresses time series andrews1993tests, bai2003computation, delgado2000nonparametric. Recent extensions to multivariate domains either impose structural assumptions park2022jump, herlands2019change, zhu2014spatially, omit regression estimation madrid2022change, require treatment status to be observed herlands2018automated, or assume known discontinuity locations cheng2023regression, narita2021algorithm, abdul2022breaking, cattaneo2022regression, keele2015geographic. Classical multivariate jump regression approaches qiu1997jump, korostelev1993minimax, o1994regularized, muller1994maximin, donoho1999wedgelets, li2017bayesian, qiu1998discontinuous often require knowing partition counts or limit dimensionality. The fused lasso literature tibshirani2005sparsity, rinaldo2009properties, harchaoui2010multiple has mostly focused on piecewise constant approximations. Recent approaches have made progress on estimating multivariate, piecewise-smooth regression surfaces with discontinuities, but do not develop any formal approaches for estimating the jump sizes and locations, instead relying on post-estimation thresholding without formal guarantees park2022jump, caroccia2020mumford, zhu2014spatially. Finally, note that there is also a rich but distinct literature on “multiple changepoints” that deals with the setting where the outcome variable is multivariate or functional aue2009break, matteson2014nonparametric, fryzlewicz2014wild, kovacs2023seeded, madrid2022change, dubey2020frechet, horvath2021monitoring, as opposed to our setting where the regressors are. A notable contribution in that literature is madrid2023change who allow the multivariate output signals to be correlated and only satisfy mixing conditions, in line with our setting for a univariate output with multivariate domain.
Our approach, on the contrary, explicitly estimates the piecewise-smooth regression function jointly with its jump locations and sizes. Moreover, we provide formal convergence and identification results for both, which to our knowledge are the first of their kind for the general multivariate setting. Additionally, all existing multivariate estimation approaches assume that the sampling noise is i.i.d. across data points or pixels, while we allow it to be correlated under a general mixing process. This is crucial for frequently occurring spatiotemporally dispersed jump processes like climate regimes, geographic policy discontinuities, or financial time series. Finally, unlike many existing methods, which often work on a fixed lattice, we allow for randomness in both the locations and the sampling error of the data points.
In this section, we develop the statistical framework for estimating regression functions with unknown discontinuities of any dimension. All mathematical notation is detailed in Appendix (ref).
We consider $n$ randomly sampled units indexed by $i=1,2,\ldots,n$ for which we observe a potentially multivariate regressor $X_i\in \mathbb{R}^d$ and an outcome of interest $Y_i\in\mathbb{R}$. The regression model is
where the $\varepsilon_i$ are unobservable error terms. Our goal is to estimate both the regression surface $f(x)$ and its discontinuity set $S_f$. To handle multivariate regressors, we need to extend the notion of a univariate discontinuity to multiple dimensions. Following Definition (ref) in Appendix (ref) ambrosio2000functions, we characterize jump points of a function $u$ on $\mathcal{X}$ through two distinct limiting values, $u^{+}(x)$ and $u^{-}(x)$, called traces. These traces reflect how $u$ behaves as one approaches a point $x$ in the domain from opposite sides, with an orientation vector $\rho$ defining these “sides." The collection of all such jump points forms the (approximate) discontinuity set $S_u$.
Standard approaches to this estimation problem face several fundamental limitations. Smoothing methods do not explicitly estimate discontinuity locations or sizes, but merely approximate the regression surface while preserving discontinuities. This leads to downward-biased jump size estimates in practice caroccia2020mumford. Simple thresholding approaches that estimate jump locations via gradient steepness almost surely misclassify discontinuities in noisy settings and fail to take into account the spatial regularity of the function that produces the jumps, resulting in severely biased estimates and high false positive and negative rates for jump detection chan2005image. The Mumford-Shah (MS) functional mumford1989optimal addresses these issues by combining smoothing and thresholding in a unified framework.
Assuming the random variable $X$ has density $f_X$, our statistical version of the MS functional takes the form:
The functional combines three components: 1) A regression term that minimizes the mean-squared error between the true function $f$ and its approximation $u$; 2) a roughness penalty that controls the gradient size of $u$ away from the discontinuity set $S_u$; 3) a penalty on the “size" of the discontinuity set, measured by its $(d-1)$-dimensional Hausdorff measure $\mathscr{H}^{d-1}$. The goal is to estimate the minimizer $u$ of this function, which immediately delivers an estimate of the discontinuity set $S_u$ as well as the jump sizes $u^+(x) - u^-(x)$ for each jump point $x \in S_u$.
The hyperparameters $\lambda, \nu\geq0$ control the balance between function smoothness and boundary complexity, with the gradient term's weight normalized to 1. The squared Euclidean norm in $\mathbb{R}^d$ is denoted by $|\cdot|^2$.
The MS functional's non-convexity through its third term threatens consistent estimation, as optimization routines need not converge to global optima. In practice, this manifests as spurious discontinuities and boundaries (brown2012completely; pock2009algorithm). A common approach to non-convex optimization is to re-estimate the model multiple times with different initializations and select the solution that minimizes the objective function. However, this approach becomes computationally infeasible for functional minimands with dimension larger than 1, as the set of potential initializations grows too large. Moreover, it is not guaranteed to converge to the true global minimum. We resolve this using calibrations alberti2003calibration, pock2010global, which lift the problem to a higher dimension by considering the graph $\Gamma_f$ of $f$. The key is finding a vector field $p:\mathbb{R}^{d+1}\to\mathbb{R}^{d+1}$ on the lifted space that maximizes flux through the graph. This yields the convex relaxation: \[ E(u)=\sup _{p \in K} \int_{\mathcal{X} \times \mathbb{R}} p \cdot D \mathds{1}_u \] with a convex constraint set
Here, $t$ indexes the lifted dimension, and $\mathds{1}_u(x,t)$ is the indicator function of the subgraph of $u(x)$, taking value 1 if $t<u(x)$ and 0 otherwise. The measure $D\mathds{1}_u\coloneqq (D_1\mathds{1}_u,\ldots, D_{d+1}\mathds{1}_u)$ is $(d+1)$-dimensional Radon measure, with $p^x$ and $p^t$ denoting the first $d$ and last components of $p\in C_0\left(\mathcal{X} \times \mathbb{R}, \mathbb{R}^{d+1}\right)$ respectively. To also achieve convexity over the target function space, we follow pock2009algorithm and optimize over the more general space of functions mapping to the unit interval:
Here, $SBV(\mathcal{X}\times\mathbb{R},[0,1])$ denotes the space of special functions of bounded variation mapping to the unit interval. This class is particularly well-suited to our problem as it naturally accommodates functions with discontinuities while still being small enough to be controllable.
This optimization problem is well-behaved due to its convexity and admits a solution under mild conditions:
While the solution $v^*$ need not take the form of an indicator function---which is required for the relaxation to be tight---we prove that it converges to one as $\lambda\to +\infty$ in a certain topology (Theorem (ref)). In practice, our automatic hyperparameter selection procedure selects values of $\lambda$ that are many orders of magnitudes larger than the selected $\nu$, and as such also leads to estimated solutions that are indicator functions in practice. Following pock2009algorithm, we obtain the underlying function $u$ by thresholding $v$ at the $0.5$-level set. When $v$ is an indicator function, any threshold $t \in (0,1)$ would yield equivalent results; the midpoint is natural when the solution deviates slightly from an indicator.
This section establishes the mathematical and statistical properties of the estimator. We first show that the solution to (ref) recovers both the true discontinuity set and jump sizes as $\lambda\to+\infty$ with fixed $\nu>0$. We then prove consistency of our empirical estimator via $\Gamma$-convergence dal2012introduction. All proofs appear in Appendix (ref).
We also provide practical tools for implementation: a data-driven method for selecting tuning parameters $\lambda,\nu$ using Stein's unbiased risk estimate (SURE), and two approaches to uncertainty quantification. The first uses subsampling with estimated convergence rates politis1999subsampling and should be used for estimation procedures, while the second employs split conformal inference lei2018distribution for computationally efficient uncertainty quantification in high-dimensional prediction.
We first establish that the true function $f$ is identified under mild regularity conditions. Specifically, we show that for fixed $\nu>0$ and as $\lambda\to+\infty$, the solution to (ref) recovers both the correct discontinuity set and jump sizes. This complements and extends previous results richardson1992limit,morini2002global which established convergence of the discontinuity set for the classical non-convex Munford-Shah functional. This result is the first to provide mathematical and statistical properties for the convexified problem.
We require two assumptions, following richardson1992limit. Let $S_f$ denote the discontinuity set of $f$, which includes the jump set $J_f$ and coincides with it $\mathscr{H}^{d-1}$-almost everywhere by the Federer-Vol'pert theorem ambrosio2000functions.
Assumption (ref) imposes standard boundedness on the functions and their jump sizes. Assumption (ref) ensures regularity both at and away from the discontinuity set. The first part implies connectedness of the discontinuity set when measured with $\mathscr{H}^{d-1}$ on $\mathcal{X}$ (e.g., precluding isolated points in two dimensions). This ensures $\mathscr{H}^{d-1}((\mathcal{X}\cap \bar{S}_f)\setminus S_f) = 0$, meaning the closure $\bar{S}_f$ coincides with $S_f$ $\mathscr{H}^{d-1}$-almost everywhere ambrosio2000functions.
The proof of (ref) requires $\lambda$ to diverge “fast enough” based on the countable rectifiability mattila1999geometry of the jump set: for any given sequence $\{r_m\}$ converging to 0 defining the lower density ambrosio2000functions of the jump set, we need $\{\lambda_m\}$ such that $\lambda_mr_m\to+\infty$. This is not a restriction since $\lambda$ can be arbitrarily large. In practice, $\lambda$ must substantially exceed $\nu$ to identify correct jump sizes. Our hyperparameter selection procedure consistently finds such values, with $\lambda$ orders of magnitude larger than $\nu$, yielding convergent estimates of both jump set and sizes as $N$ grows. We now turn to the practical implementation and the corresponding estimator.
We now develop the empirical version of the estimator for observed data. Consider a random sample $\{X_i,Y_i\}_{i=1}^n$ from distribution $P_{Y,X}$, with errors $\varepsilon_i$, which can be correlated and not identically distributed, as long as they satisfy a law of large numbers. This is important for the spatial setting in our application, where the error terms in the nonparametric regression are correlated in general. Sufficient for this are mixing conditions. In the following, to guarantee consistency of the kernel density estimator $\hat{f}_{X,h(n)}$, we assume that the marginal data generating process $\{X_i\}_{i\in\mathbb{N}}$ is stationary. This can be relaxed irle1997consistency.
The key challenges in moving to random data are twofold:
To accommodate the random locations, we discretize (ref) on a regular grid following pock2009algorithm. Specifically, we overlay an $(N_1 \times N_2 \times \ldots \times N_d) \times S$ pixel grid on the random point cloud $\{X_i\}_{i=1}^n$, defining $ \mathcal{Q}_N=\{(k_1,\ldots, k_d, s): k_i=1,2, \ldots, N_i; \ldots; k_d = 1,\ldots N_d; s=1,2, \ldots, S\}. $ To ease notation, we assume without loss of generality that all dimensions are discretized in the same manner, $s\coloneqq k_{d+1}$, $N_1,\ldots,N_{d+1} = N$. Also, denote $k\coloneqq (k_1,\ldots,k_{d+1})$, and $k_{-}\coloneqq (k_1,\ldots,k_d)$. The empirical analogue of $v\in C$ becomes \[\hat{v}_{Nn}(x)\equiv \hat{v}_{Nn}(\tilde{x},t)\coloneqq \sum_{1\leq k_1,\ldots,k_{d+1}\leq N}v_{k_1,\ldots,k_{d+1}}\mathds{1}\left\{x\in Q_{k_1,\ldots,k_{d+1}}\right\}\] where $v_k$ is assigned to the center of cube $Q_k$. We define the empirical analogue $p_N$ of $p\in K$ to take values on the boundary of the respective cubes. That is, for each $j=1,\ldots, d+1$, the corresponding value $p_{k_1,\ldots,k_{j}+\frac{1}{2},\ldots,k_{d+1}}$ lies on the boundary $\partial Q_{k_1,\ldots,k_{j},\ldots,k_{d+1}}\cap \partial Q_{k_1,\ldots,k_{j}+1,\ldots,k_{d+1}}$. This preserves duality in the empirical setting, as the empirical analogues of $p$ are defined on boundaries of the cubes, while the analogues of $v$ are defined in the center chambolle2021learning. Problem (ref) becomes
with
where $D_Nv$ is the normalized forward difference and $\langle p, Dv\rangle_N$ the scalar product between two vectors, \[\langle p_N,D_Nv_N\rangle_N =\sum_{0\leq k_1,\ldots, k_{d+1}\leq N} N\left(v^\uparrow_{k_1,\ldots, k_{d+1}} - v_{k_1,\ldots, k_{d+1}}\right) p_{k_1,\ldots,k_{d+1}}, \] with $v^\uparrow_{k_1,\ldots,k_j,\ldots, k_{d+1}}\coloneqq v_{k_1,\ldots,k_{j}+1,\ldots,k_{d+1}}$ and $p_k^\uparrow = p_k\cdot e^\uparrow$. The constraints hold for all grid points $k$, and $s_1 \leq k_{d+1} \leq s_2$ iterates over indices with $1\leq s_1 \leq k_{d+1} \leq s_2 \leq S$.
For function estimation, we use weighted averages within each grid cube: \[\hat{f}_{Nn}(\bar{x}_{k_{-}})=\sum_{i: X_i\in Q_{k_{-}}} w_i Y_i\] where $\bar{x}_{k_{-}}$ is the cube center and weights $w_i$ sum to 1; for simplicity, we assume uniform weights, i.e., $w_i = 1 / \#\left\{i: X_i\in Q_k\right\}$. We estimate the density $f_x$ via a standard Nadaraya-Watson estimator \[\hat{f}_{X,h(n)} = \frac{1}{nh(n)^d}\sum_{i=1}^n K\left(\frac{X_i-x}{h(n)}\right),\] where $h(n)$ is a bandwidth and $K(\cdot)$ is a multivariate kernel. We assume that $K:\mathbb{R}^d\to\mathbb{R}$ satisfies $K(x\rho) \geq K(x)$ for all $x\in\mathbb{R}^d$ and $0\leq\rho\leq 1$, $\int_{\mathbb{R}^d} K(x) dx = 1$, $K(x) = O(|x|^{-n-\rho})$ as $|x|\to+\infty$ for some $\rho>0$, and that $\sup_{x\in\mathbb{R}^d} K(x) <+\infty$. In the following, we call such a kernel regular irle1997consistency.
We make the following assumption on the data-generating process.
These regularity assumptions are standard. Part (i) requires a bound on the mixing parameter $s$ that is very minor, since $r>1$ can be as close to unity as we want. It involves the rate $m$ at which the number $N^d$ of grid points grows with the number of data points $n$. The slower this growth, the smaller the mixing constant can be. The strict stationarity assumption on $\{X_i\}_{i=1}^n$ is done for convenience in the proof and can be relaxed. Further, we allow the penalty term $\lambda_n$ to vary with the data, but very slowly, in order to still obtain the almost sure convergence. The moment conditions in Part (ii) are very minor and allow for unbounded error terms $\varepsilon_i$ as long as their tails decay quickly enough. The regularity for kernel and histogram estimators in Part (iii) are standard and taken from the kernel estimator results in irle1997consistency.
Before proving the full statistical result, we provide a consistency result for the simpler setting where we pretend $f(x)$ is known and does not have to be estimated. This result is novel in the image recognition literature and complements recent convergence results in this area chambolle2021learning. It also showcases the complication when introducing the two sources of randomness.
Now we add the randomness to $f$ as well as $f_X$. The fact that $f$ has to be estimated through random draws complicates the estimation procedure significantly. In fact, any estimator $\hat{f}_{Nn}$ for $f$ that averages out the error caused by $\varepsilon_i$, including the proposed histogram estimator, converges to a precise representative $\tilde{f}$ of $f$. This means that it is consistent only for points $x\not\in S_f$. For points $\tilde{x}\in S_f$, $\hat{f}_{Nn}(\tilde{x})$ converges almost surely to $\tilde{f}(\tilde{x}) = (1-\theta) f^+(\tilde{x}) +\theta f^{-}(\tilde{x})$, that is, a convex combination of the values of the traces $f^+(\tilde{x})$ and $f^-(\tilde{x})$. The weight $\theta$ depends on the angle of the hyperplane induced by the approximate differential and the weighting from $f_X$. If all $Q_k$ are cubes and $f_X$ is the density of the uniform distribution, then $\theta=1/2$. Importantly, $0<\theta<1$, which implies that the proposed estimator identifies the location of the jumps of $f$ perfectly. This is captured in the following result.
Complementing this with a compactness result yields convergence of the minimizer:
This shows that the estimator is consistent for a precise representative of $f$, identifying the correct jump locations as $\lambda_n\to+\infty$. As shown in the simulations in Figure (ref), the proposed estimator also allows to identify the correct jump size, not just the location. The key for this is the weak$^*$-convergence of the optimizer as proved in Theorem (ref): by the closure of SBV ambrosio2000functions, this implies that the traces of $v_n$ converge $\mathcal{H}^{d-1}$-almost everywhere to the traces of $v$. Since the traces are identified, this also identifies the jump sizes directly. Finally, standard convergence results for level-sets camilli1999note imply that the $0.5$-level set, which is our estimator for the true function, also converges almost surely.
In practice, following the vision literature strekalovskiy2014real, we detect the empirical discontinuity set by thresholding the estimated gradient, \[ \hat S_u \;=\; \bigl\{x\in\mathcal X \;:\; \|\nabla\hat u_{Nn}(x)\|\;\ge \sqrt{\nu}\bigr\}. \] That gives us the estimated jump locations. To estimate the corresponding jump sizes, for every boundary voxel \(x\in\hat S_u\), we march in the estimated outward normal direction \(\hat{\rho}_{N_n}(x)\) and in the inward direction \(-\hat{\rho}_{N_n}(x)\) until we reach the first voxels \(x_{N_n}^{+}(x)\) and \(x_{N_n}^{-}(x)\) that lie outside the boundary mask \(\hat S_u\). The empirical jump–size estimator is \[ \widehat{\Delta}_{N_n}(x) \;=\; \hat u_{N_n}\!\bigl(x_{N_n}^{+}(x)\bigr) \;-\; \hat u_{N_n}\!\bigl(x_{N_n}^{-}(x)\bigr), \] i.e.\ the difference between the “pure” function values immediately to the right and left of the interface. Because the march follows \(\hat{\rho}_{N_n}(x)\), the construction automatically captures the correct orientation, so no grid refinement is needed. Then, Corollary (ref) implies convergence of the empirical normal $\hat{\rho}_{N_n}(x)$ to the population normal to the jump set, and hence of the estimated jump sizes to their population counterparts.
Finally, we construct confidence bands using subsampling politis1999subsampling, which we choose over bootstrap methods due to the estimator's non-smooth nature. The Online Appendix provides implementation details and describes a more computationally efficient alternative using conformal prediction.
While Theorem (ref) establishes that the population estimator recovers the true jump locations and sizes as $\lambda$ grows for fixed $\nu$, optimal hyperparameter choice in finite samples depends on the data. We select $\lambda,\nu$ by minimizing Stein's unbiased risk estimate (SURE) stein1981estimation, which provides an asymptotically unbiased estimate of mean-squared error under Gaussian errors.
Given estimator $\hat{u}_\theta(Y)$ with hyperparameters $\theta = (\lambda, \nu)$, the Stein estimator is
where $\operatorname{div}_Y \hat{u}_\theta(Y) \coloneqq \sum_{i=1}^N \frac{\partial \mathrm{\hat{u}_{\theta,i}}(Y_i)}{\partial Y_i}$ is the divergence with respect to the data. For continuous and bounded operators with well-defined divergence, $\eta(\hat{u}_\theta)$ is unbiased for $MSE(\hat{u}_\theta(Y)) \coloneqq \frac1N \lvert \hat{u}_\theta(Y) - f \rvert^2$ ramani2008monte.
Since the estimator lacks a closed-form divergence, we use a Monte Carlo approximation by perturbing the data with random noise $\mathbf{b} \sim N(0, 1)$:
Following lucas2022hyperparameter, we average over $R$ draws of $\mathbf{b}^{(r)}$ to obtain $\bar{\eta}^R(\hat{u}_\theta(Y)) \coloneqq \frac1R \sum_{r=1}^R \eta_{\delta, \mathbf{b^{(r)}}}(\hat{u}_\theta(Y))$, setting $\delta=0.01$. Monte Carlo simulations in the next section confirm that this approach selects effective hyperparameters.
{
}
We validate the estimator on 1D, 2D, and 3D simulations with piecewise-smooth functions featuring jumps of varying magnitudes (Cohen’s \(d = 0.25,\) \( 0.5,\) \(0.75\)). Additive Gaussian noise is set to \(\sigma = 0.05\). We discretize the lifted dimension into 32 points pock2009algorithm and select hyperparameters \(\lambda,\nu\) via a Monte Carlo version of Stein’s unbiased risk estimate (MC-SURE) ramani2008monte, lucas2022hyperparameter on a \(20 \times 20\) grid.
Figure (ref) shows the true functions, noisy point clouds, and estimated surfaces for sample sizes of \(n = 5{,}000\) (1D), \(n = 10{,}000\) (2D), and \(n = 50{,}000\) (3D). The method accurately recovers both location and magnitude of discontinuities. In the 1D case, we also plot 95% confidence bands computed via subsampling politis1999subsampling, showing that all estimated jumps are significant.
Table (ref) reports mean squared error (MSE), jump-size bias, and misclassification rates for boundary detection, averaged over 300 (1D) and 100 (2D) simulations. All measures improve as \(n\) grows, confirming our convergence results -- except for a slight and mechanical increase in the false negative rate for the jump locations due to the finer grid sizes. The false positive and negative rates for the estimated jump locations rapidly decrease to zero, indicating highly accurate jump location detection even in small samples. The estimated jump sizes for the 2D case converge to estimates that different by less than 0.5% from the true jump sizes, except for the smallest jumps, which still converge but at a slightly slower rate. Notably, the data-driven selection of hyperparameters consistently finds \(\lambda \gg \nu\), aligning with Theorem (ref).
Internet shutdowns—deliberate disruptions of internet or electronic communications—have recently drawn global attention, reaching a record high in 2022 rosson2023weapons. India, with its burgeoning digital economy, has implemented at least 646 such shutdowns between 2018 and 2023 InternetShutdowns2024, more than any other country, primarily to quell protests, communal violence, and cheating in examinations hrw. In this section, we exploit the geographic discontinuities caused by these shutdowns to estimate their short-term economic impact.
The FDR estimator is well-suited for this purpose because the areas exposed to the shutdown—where internet connectivity abruptly dropped—were unknown a priori. Its multidimensional formulation allows us to explicitly estimate the geography of the targeted areas. This is of independent interest, to understand what areas were affected, but also enables us to precisely estimate the impact of the shutdown by looking at the size of the jumps at the estimated area boundaries.
We focus on a shutdown imposed by the Rajasthan state government on September 26, 2021, to prevent cheating on the Rajasthan Eligibility Exam for Teachers. This high-stakes exam, not held since 2018, drew hundreds of thousands of candidates in 2021 amid heightened concerns of cheating after several scandals Purohit. In the days before the exam, multiple district governments announced a mobile internet shutdown from 6 am to 6 pm.\footnote{Specifically, the District Magistrates or Divisional Commissioners of Ajmer, Jhunjhunu, Kota, Bundi, Baran, Jhalawar, and parts of Udaipur issued suspension orders Mishra_2021.}
Several other districts also cut telecom services without notice, ultimately affecting an estimated 25 million people yeung2021internet. Because the Rajasthan telecom circle does not perfectly align with administrative boundaries, and connectivity depends on how cell towers, service providers, and mobile devices interact, it was unclear which districts—and what parts of those districts—would be affected. Where the shutdown was enforced, it was complete, as local governments compelled providers to disable both transmission (cell towers) and point-of-contact (mobile phones). This meant users could not connect to towers outside the shutdown area, inducing discontinuous drops in connectivity.
Although the shutdown technically targeted mobile services, it effectively amounted to a near-complete internet blackout. In 2022, Wi-Fi accounted for just 0.08% of wireless connectivity in India, and only 3.74% of internet subscribers had wired connections trai2023indian. This disruption reverberated through the economy. In Jaipur alone, 80,000 shops reportedly closed economist_shutdowns, due to the ubiquity of mobile-based payment systems like UPI or Google Pay. Although cash is accepted, it is less convenient in a system where even small transactions often rely on mobile platforms nyt_mobile_payments, leading many ATMs to run dry during similar shutdowns jaipur_shutdown. A large share of retail, hospitality, and mobile-app-based transactions also ground to a halt.
More generally, these shutdowns directly disrupt India’s digital economy, which accounted for about 22% of GDP in 2019 and continues to expand. While wired connections remained available, losing mobile connectivity still undermined supply-chain tracking, process automation, distribution networks, remote work, and customer support—especially for small and rural businesses that rely entirely on mobile internet kathuria2018anatomy.
A commonly cited framework for estimating the costs of internet shutdowns is the back-of-the-envelope method from west2016internet,\footnote{See also top10vpn2023 and \url{https://netblocks.org/projects/cost}.} which uses data on the digital economy’s size, mobile penetration, and the “digital multiplier” quelch2009quantifying. While transparent, this approach does not measure actual local economic impacts, and there appear to be no academic estimates of these, despite shutdowns’ global prominence. In development economics, however, research on digital infrastructure expansions suggests notable economic benefits. roller2001telecommunications, for instance, found that each percentage-point rise in mobile penetration raises economic growth by around 0.15%. Other studies show welfare improvements, lower living costs, and increased employment from mobile rollouts jensen2007digital, bjorkegren2022network, couture2021connecting, zuo2021wired, hjort2019arrival, consistent with digitization reducing search, replication, transportation, tracking, and verification costs goldfarb2019digital.
In contrast, an internet shutdown disrupts those welfare-improving channels, creating economy-wide ripple effects. Moreover, losing established connectivity is not simply the inverse of gaining it; digital infrastructures, once in place, play “an ‘enabling function’ across all critical infrastructure sectors” CISA. Our findings therefore complement, rather than mirror, existing work on connectivity expansions.
Lastly, using mobile device data to proxy socioeconomic activity is well-established in fields such as remote sensing, network science, and complex systems vscepanovic2015mobile, kung2014exploring, frias2012relationship, eagle2010network, blumenstock2015predicting, as well as in industry naef2014using and government worldbank2022ukraine. More recently, economists and policy researchers have adopted these data to study social and economic questions kreindler2021measuring, van2023public.
We analyze anonymized device-level location data from data provider Veraset to assess the effects of the shutdown at a fine spatio-temporal granularity. The data consist of “pings”, timestamped GPS locations shared by mobile devices with apps that use common software development kits. Inside our sample area around Rajasthan, we observe 126 million unique pings from 3.8 million devices during the four Sundays of September 2021 between 6 am and 6 pm, the shutdown time window. To capture the disruption to mobile connectivity, we define
where $i$ indexes a 5$\times$5km grid cell overlaid on Rajasthan; $t_0$ indexes September 26, 2021, between 6 am and 6 pm (the shutdown period); and $t_0-3,\ldots,t_0-1$ index the same time window on the three previous Sundays. $Pings_{it}$ counts mobile device pings in grid cell $i$ in period $t$. The normalization by the prior month's Sunday average removes regional differences in absolute activity levels.
To estimate economic effects, we compute an analogous measure $\overline{Econ}_i$ that captures activity around commercial points of interest (POIs). We combine two datasets --- SafeGraph Places (global retail chains) and OpenStreetMap (local businesses) --- to identify approximately 108,000 retail, commercial, and industrial locations (see Appendix (ref) for detailed classifications). We classify a ping as economically relevant if it originates within 150 meters of a POI. The economic measure is cast onto a coarser 40$\times$40km grid intersected with municipalities to ensure sufficient observations per cell.
Some devices continue emitting signals during the shutdown because location data is determined via satellite connections, which were unaffected. Apps with background location permissions cache these positions for later transmission, enabling measurement of economic activity even in shutdown areas. The representativeness of this cached location data and additional robustness checks are discussed in Appendix (ref).
We apply the FDR estimator to the data described above, starting with the mobile device signal in Figure (ref).
\paragraph*{Mobile signal} The raw data (Figure (ref)) shows a stark drop-off in signal that aligns with the Rajasthan state boundary, with an unexpected pocket of activity in the northwest intersecting the districts of Jaisalmer, Jodhpur, and Nagaur. Areas without prior signal are normalized to 100, while areas unaffected by the shutdown show signals averaging 130% of the Sunday September average. This reflects both mechanical noise from low spatial granularity and the broader economic recovery in India at that particular time in the COVID-19 pandemic Bhaduri_2021.
For hyperparameter selection, we search over $\lambda \in [1,100]$ and $\nu \in [0.001, 0.1]$ to account for the large amount of noise in the mobile data. We winsorize at the 90th percentile to handle outliers from temporary gatherings and set $f_X$ to the uniform distribution after verifying robustness to alternatives. The resulting regression function $\hat{u}$ (Figure (ref)) is smooth within the two partitions, with the estimated jump set $S_u$ precisely capturing the state boundary and the unexpected area of continued connectivity within the state. The effect curve ((ref)) shows an average reduction in mobile signal of 100% relative to baseline, with higher reductions near Gujarat and Madhya Pradesh. About 25% of the baseline signal remains inside the shutdown area due to satellite-based location caching, enabling our subsequent economic analysis. Overall, we estimate that the shutdown was highly effective in its goal of disrupting connectivity, inducing a large drop in the mobile device signal throughout the targeted area.
\paragraph*{Economic activity} Figure (ref) shows the results for economic activity. We construct 95% confidence bands around the function gradient using 200 subsamples at 4 different subsample sizes, estimating a convergence rate of approximately $\sqrt{n}$. Hyperparameters are selected from $\lambda \in [1,5]$ and $\nu \in [0.075, 0.15]$, adjusted for the coarser grid. Despite lower resolution, the estimated jump set closely tracks the mobile signal discontinuities, overlapping with state boundaries and identifying the northwestern pocket of connectivity.
The shutdown reduced economic activity by approximately 35% in the eastern region and 60% in the western region, though the western estimate is less reliable due to sparse data and proximity to Pakistan. Outside the shutdown area, activity increased 0-5%, consistent with India's broader economic recovery at the time woloszko2020tracking. Scaling our estimate by 0.7—the most conservative estimate of mobility-GDP correlation from the literature dong2017measuring, frias2012relationship, spelta2021mobility—yields a 25% reduction in economic activity. This substantially exceeds previous estimates of shutdown costs, which predicted impacts around 16.4% of annualized GDP.
This discrepancy likely reflects two factors: previous estimates focused narrowly on digitized economic activity, missing spillovers to non-digital sectors, and relied on outdated digital economy multipliers from 2009 quelch2009quantifying that understate digital infrastructure's current importance, particularly in rapidly digitalizing economies like India.
\paragraph*{Implications} Our results reveal that arbitrary internet shutdowns inflict larger short-term economic damage than previously understood. While some recovery might occur—though we find no evidence of this in mobile or economic activity the following day (Figure (ref))—the magnitude of our estimates underscores digital infrastructure's critical role in modern economies. The findings highlight an important asymmetry: while the benefits of expanding internet access accrue gradually, the costs of removing it are immediate and substantial once digital infrastructure becomes integrated into economic activity. This emphasizes not just the importance of expanding internet access, but of ensuring its continued reliability UN2022.
{
We introduced Free Discontinuity Regression (FDR): the first fully non-parametric framework that simultaneously estimates a multivariate regression surface and the geometry—and magnitude—of its unknown discontinuities. Built on a statistical reinterpretation of the convexified Mumford–Shah functional, FDR overcomes three long-standing bottlenecks in high-dimensional jump surface estimation:
A data-driven version of Stein's Unbiased Risk Estimate selects the two tuning parameters automatically and in a way that is consistent with our identification theorems. Large-scale simulations from 1D to 3D confirm the predicted convergence of both surfaces and jump sizes, while an application to internet shutdowns in India showcases FDR’s empirical power—revealing a 25–35% drop in economic activity, far larger than headline estimates. This points to a large asymmetry in the effects of internet expansions and shutdowns, underscoring the internet's critical role in modern economies.
\paragraph*{Outlook} By combining theoretical guarantees with efficient optimization procedures, the proposed estimator provides a robust and practical statistical estimation tool for multivariate changepoint estimation. Immediate extensions include (i) dynamic or panel settings where jump surfaces evolve over time, through vectorial coupling of the multiple output channels strekalovskiy2012convex, (ii) adaptive or multiresolution grids for ultra-high-resolution imagery, (iii) instrumental-variable and causal designs that leverage the estimated discontinuity set as an endogenous boundary, and (iv) computational improvements through statistical versions of the sublabel-accurate approaches in the image literature mollenhoff2017sublabel. The hope is that the framework will spur new research into the estimation and application of jumps in high-dimensional processes and further bridge the gap between computer vision and statistical estimation.