Spectral Bounds and Exit Times for a Stochastic Model of Corruption

1. Introduction

Institutional corruption, understood as the abuse of public power for private benefit, represents a critical obstacle to socioeconomic development and a threat to democratic governance. In recent years, it has been recognized that stochastic mathematical models—particularly those inspired by epidemic dynamics (SIR-type models)—offer a promising framework to analyze the spread of corruption and to evaluate the impact of anti-corruption policies; see, for example, [1,2,3,4]. These models often incorporate random factors (Brownian or Lévy noise) that capture intrinsic uncertainties, such as the random detection of offenses or variability in the effectiveness of sanctions. The stochastic formulation allows, among other things, the quantification of exit times (i.e., the first hitting times of acceptable corruption levels) via Kolmogorov/Fokker–Planck type partial differential equations, which are essential to assess the temporal effectiveness of governmental interventions.

In this work, we introduce a novel stochastic model for institutional corruption aimed at evaluating public policies. We start from a framework that explicitly includes policy-related variables such as investment in transparency, ethical education, and legal sanctions. By incorporating Brownian noise into the transition rates, we capture unpredictable variability in the effectiveness of such measures and in the spread of corrupt behavior. Mathematically, this leads to a system of multivariate stochastic differential equations (SDEs). A distinctive feature of our model is the emphasis on exit times of corruption under different policy scenarios. We define a domain of interest D (for instance, high levels of corruption) and compute the first exit time ${\tau }_D$ to more acceptable levels. The distribution and density of ${\tau }_D$ are obtained analytically from the Kolmogorov/Fokker–Planck equation, whose numerical solution is computed using FreeFEM. This enables, for example, the evaluation of how changes in audit funding or educational campaigns affect the expected time to achieve a “reduction” of corruption.

The proposed formulation addresses several issues identified in the literature, see [14]. First, it explicitly integrates public policy factors into the stochastic dynamics, thereby connecting model parameters to real interventions. Second, by focusing on exit times, it provides a clear temporal metric (rather than purely population-based measures) to assess policy efficiency. Finally, it employs robust numerical methods to solve the PDEs associated with the exit time ${\tau }_D$.

A central theoretical result of this work arises from the study of the linearization of the system around the deterministic equilibrium. We show that bounds on the second-order moment (the mean square of perturbations) are governed by the spectral structure of the linearized matrix—that is, by its eigenvalues rather than by the magnitude of a single noise term. In particular, the presence of a dominant positive eigenvalue leads to exponential growth of variability and limits the temporal validity of any linear approximation, rendering its conclusions essentially local in time. Conceptually, this explains why small variations in institutional reactivity parameters or in the effectiveness of sanctions can amplify uncertainty in non-intuitive ways, producing “overcorrection” effects. From a public policy perspective, this suggests that interventions must be designed with the collective interaction of parameters in mind (i.e., the system’s spectral structure), rather than focusing solely on reducing volatility in a single component.

We apply the proposed framework to the case of Mexico and compute exit times numerically for different policy calibrations. The results show that persistence in undesirable states does not depend on a single parameter in isolation but rather on the interaction between institutional laxity, punishment effectiveness, and governmental reactivity: certain combinations can significantly prolong high levels of corruption. Moreover, the volatility of institutional laxity affects exit speed in a non-monotone way: reducing policy uncertainty often accelerates recovery, but excessively reactive responses without effective sanctions may induce overcorrections and increase uncertainty. Based on our analysis for Mexico, we recommend prioritizing interventions that modify the system’s drift (e.g., strengthening the effectiveness of sanctions and implementing institutional reforms that reduce regulatory laxity), rather than relying on random fluctuations; in addition, institutional responses should be graduated to avoid overcorrections that exacerbate uncertainty. It is essential to incorporate time-based metrics derived from exit times to establish realistic policy horizons and evaluation criteria, thereby identifying high-risk parameter combinations. Such measures allow for faster and more robust improvements in reducing corruption while minimizing unpredictable effects stemming from complex interactions among laxity, sanction, and reactivity.

The organization of the paper is as follows. In Section 2, we present the deterministic model introduced in [14], which we then perturb randomly. Section 3 introduces the stochastic dynamics and establishes global existence. Section 4 studies a linear approximation of the stochastic system. Section 5 and Section 6 analyze the exit time through partial differential equations and their numerical solution. Finally, in Section 7 we apply the model to the case of corruption in Mexico.

2. A Deterministic Model for Corruption

To model corruption we consider the following system of ordinary differential equations

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dC\left(t\right)}{dt}}& =& \alpha C\left(t\right)P\left(t\right)-{\gamma }_{1}C\left(t\right)L\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dL\left(t\right)}{dt}}& =& -\beta L\left(t\right)+{\gamma }_{2}C\left(t\right)L\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dP\left(t\right)}{dt}}& =& \delta P\left(t\right)(\xi -P(t\left)\right)-ϵC\left(t\right),\hfill \end{array}$$

where the variables represent:

• $C\left(t\right)$: Normalized corruption perception index ($C=0$ no corruption, $C=1$ total corruption).
• $L\left(t\right)$: Proportion of corrupt acts that have been prosecuted and sanctioned by law ($L=0$ total impunity, $L=1$ full justice, zero impunity).
• $P\left(t\right)$: Laxity of anti-corruption policies ($P=0$ strict and effective, $P=1$ lax or non-existent).

All parameters of the model are positive real numbers and have the following interpretation:

• $\alpha$ (corruption growth induced by institutional laxity): Represents the rate at which corruption perception grows when corruption already exists and anti-corruption policies are permissive. The larger $\alpha$ is, the stronger the negative impact of institutional laxity on the social perception of corruption.
• $\beta$ (institutional deterioration): Represents the natural tendency of judicial institutions to lose effectiveness in the absence of external stimuli. It is the rate at which the legal system loses its ability to punish, even without changes in corruption.
• ${\gamma }_1$ (corrective effect of punishment): Measures the ability of the judicial system to reduce the perception of corruption through punishment. A large ${\gamma }_1$ means punishments effectively reduce the social perception of corruption.
• ${\gamma }_2$ (public pressure to punish): Models how high corruption perception triggers legal action. The larger ${\gamma }_2$, the more sensitive the judicial system is to public perception, and the higher the proportion of sanctioned corrupt acts.
• $\delta$ (natural growth of laxity): The rate at which the laxity of anti-corruption policies tends to increase, modeled by logistic dynamics. It represents institutional deterioration in the absence of social or political pressure.
• $\xi$ (maximum capacity of permissiveness): The value towards which $P\left(t\right)$ tends in the absence of corrective mechanisms. It represents the maximum structurally tolerated level of laxity in anti-corruption policies.
• $ϵ$ (governmental response to corruption perception): Measures the sensitivity of policy tightening in response to corruption perception. The larger $ϵ$, the faster and more significant the institutional response when perceived corruption increases.

This deterministic model, considered in [14], captures complex dynamics among social perception, judicial efficiency, and institutional response. The parameters can be calibrated to evaluate different public policy scenarios and their effects on perceived corruption. In Section 7 we use the stochastic version of this model to study corruption in Mexico.

It is well known that a randomly perturbed model can behave quite differently from the original. However, if our premise is that the deterministic model reflects reality in some sense, then it is useful, so we want our stochastic results not to deviate too much from it. To make this more precise, note that the equilibrium points of the deterministic system are

$$(0,0,0),(0,0,\xi )\mathrm{and}{E}^{*}=(C^{*},L^{*},P^{*}),$$

where

$$C^{*}={\displaystyle \frac{\beta }{{\gamma }_2}},L^{*}={\displaystyle \frac{\alpha }{{\gamma }_1}}·P^{*},$$

(1)

and $P^{*}$ is a solution of the quadratic equation

$$\delta P^2-\delta \xi P+{\displaystyle \frac{ϵ\beta }{{\gamma }_2}}=0.$$

(2)

It is important to point out that among the four equilibrium points, the only asymptotically stable equilibrium is $E^{*}=(C^{*},L^{*},P^{*})$, where

$$P^{*}={\displaystyle \frac{\xi -\sqrt{{\xi }^2-4·\frac{ϵ\beta }{\delta {\gamma }_2}}}{2}},$$

(3)

whenever

$$\xi >2\sqrt{{\displaystyle \frac{ϵ\beta }{\delta {\gamma }_2}}}.$$

(4)

That is, $P^{*}$ is the negative root of (2).

Taking this result from the deterministic model, namely that $E^{*}$ is the only asymptotically stable equilibrium point, in what follows we are interested in studying the behavior of the stochastic model, in the next sections, around $E^{*}$.

3. A Stochastic Corruption Model

In the deterministic corruption model of Section 2 we assume that there is uncertainty in the value of the following parameters: corrective effect of punishment, public pressure to punish, and governmental response to corruption perception. Thus, if we denote by ${\gamma }_1$, ${\gamma }_2$, and $ϵ$ the corresponding average values of these parameters, then the new parameter values can be expressed as

$${\gamma }_1+{\sigma }_1{\xi }_1\left(t\right),{\gamma }_2+{\sigma }_2{\xi }_2\left(t\right)\mathrm{and}ϵ+{\sigma }_3{\xi }_3\left(t\right),$$

where ${\xi }_1\left(t\right)$ denotes white noise and ${\sigma }_i>0$ measures the dispersion of such noise, see [13]. Denoting symbolically ${\xi }_{i}dt$ by $d{W}_i\left(t\right)$, the deterministic system of Section 2 can be written in the form

$$dX\left(t\right)=f\left(X\right(t\left)\right)dt+\sigma \left(X\right(t\left)\right)dW\left(t\right),$$

(5)

where $X\left(t\right)={(C\left(t\right),L\left(t\right),P\left(t\right))}^{\top }$. Moreover, the vector-valued drift function f is

$$f=\left(\begin{array}{c}\hfill f_1\\ \hfill f_2\\ \hfill f_3\end{array}\right):=\left(\begin{array}{c}\hfill \alpha CP-{\gamma }_{1}CL\\ \hfill -\beta L+{\gamma }_{2}CL\\ \hfill \delta P(\xi -P)-ϵC\end{array}\right)$$

(6)

and the diffusion matrix $\sigma$ is given by

$$\sigma =\left({\sigma }_{ij}\right):=\left(\begin{array}{ccc}\hfill -{\sigma }_{1}CL& 0& 0\hfill \\ \hfill 0& {\sigma }_{2}CL& 0\hfill \\ \hfill 0& 0& -{\sigma }_{3}C\hfill \end{array}\right).$$

Equation (5) is a stochastic differential equation, to be understood in the Itô sense. Furthermore, $W\left(t\right)={(W_1\left(t\right),W_2\left(t\right),W_3\left(t\right))}^{\top }$ is a Brownian motion in ${\mathbb{R}}^3$ defined on a complete filtered probability space $(\Omega ,\mathcal{F},\left\{{\mathcal{F}}_t\right\},\mathbb{P})$. By $\mathbb{E}$ we denote expectation with respect to the probability measure $\mathbb{P}$.

Our first goal is to prove that the stochastic differential equation (5) has a unique global solution. Since the quantities involved in X are averages, the region of interest will be the set

$$D=\left\{{(C,L,P)}^{\top }\in {\mathbb{R}}_{+}^3:0\le C\le 1,0\le L\le 1,0\le P\le 1\right\}.$$

Theorem 1.

For any initial value ${(C\left(0\right),L\left(0\right),P\left(0\right))}^{\top }\in D$, equation (5) has a unique positive solution $X\left(t\right)={(C\left(t\right),L\left(t\right),P\left(t\right))}^{\top }\in D$, almost surely for all $t\ge 0$.

Proof. 

Since the functions $f_i$ and ${\sigma }_{ij}$ are polynomial functions, they satisfy the local Lipschitz condition and, by the theory of Itô SDEs, see [12], there exists a unique local solution $X\left(t\right)$ up to an explosion time ${\tau }_e>0$. For each integer $k\ge 1$ we define the stopping time

$${\tau }_k=inf\left\{t\in [0,{\tau }_e):min\{C\left(t\right),L\left(t\right),P\left(t\right)\}\le 1/k\mathrm{or}max\{C\left(t\right),L\left(t\right),P\left(t\right)\}\ge k\right\}.$$

Then ${\tau }_k\uparrow {\tau }_e$ as $k\uparrow \infty$. Our aim is to prove that ${\tau }_e=\infty$ almost surely, which will imply the global existence of X and that $X\left(t\right)\in D$ almost surely for all $t\ge 0$.

Consider the function

$$V(C,L,P)=(C-1-lnC)+(L-1-lnL).$$

Applying Itô’s formula to $V\left(C\right(t),L(t),P(t\left)\right)$ yields

$$dV=\mathcal{L}Vdt-{\sigma }_1(C-1)Ld{W}_1\left(t\right)+{\sigma }_2(L-1)Cd{W}_2\left(t\right),$$

where

$$\mathcal{L}V=(C-1)(\alpha P-{\gamma }_{1}L)+(L-1)(-\beta +{\gamma }_{2}C)+{\displaystyle \frac{1}{2}}({\sigma }_1^2{L}^2+{\sigma }_2^2{C}^2).$$

Since the function $\mathcal{L}V$ is continuous in D, there exists a constant $K>0$ such that

$$\mathcal{L}V\le K,\mathrm{for} \mathrm{all}{(C,L,P)}^{\top }\in D.$$

Integrating $dV$ from 0 to ${\tau }_k\wedge T$ and taking expectation gives

$$\mathbb{E}\left[V\left(X({\tau }_k\wedge T)\right)\right]=V\left(X\left(0\right)\right)+\mathbb{E}\left[{\int }_0^{{\tau }_k\wedge T}\mathcal{L}V\left(X\left(s\right)\right)ds\right]\le V\left(X\left(0\right)\right)+KT.$$

If there existed an $\epsilon >0$ such that $\mathbb{P}[{\tau }_k\le T]\ge \epsilon$, then on that event one of the components would equal $1/k$ or k, forcing $V\left(X\left({\tau }_k\right)\right)\ge lnk-1$. Thus

$$V\left(X\right(0\left)\right)+KT\ge \epsilon (lnk-1),$$

and hence

$$lnk\le {\displaystyle \frac{V\left(X\right(0\left)\right)+KT}{\epsilon }}+1,$$

which is impossible as $k\uparrow \infty$. This contradiction shows that ${\tau }_e=\infty$ almost surely. Therefore, the solution $X\left(t\right)$ is global, positive, and remains in the bounded region $\{0\le C,L,P\le 1\}$, as desired.    □

4. A Semilinear Stochastic Model for Corruption

Recall that a desirable aspect is the study of the stability of system (5). However, there is a problem: the only critical point of the stochastic system is ${(0,0,0)}^{\top }$, which is of little interest due to the utopian nature of such a political system. Moreover, it is worth noting that it is common for perturbed models (see [11]) to exhibit this type of behavior, namely, that the critical points of the deterministic system are not necessarily the same as those of the stochastic system, and when the system is randomly perturbed it either becomes extinct or explodes in finite time (i.e., it leaves the region of interest).

Next, we will study a first-order approximation of the stochastic system (5) and we will see that noise, however small, immediately implies that $E^{*}$ is no longer an equilibrium point of the system. Furthermore, we will show that the mean square of the linearized system exhibits exponential growth, in stark contrast with the behavior of the process $X\left(t\right)$ and of the deterministic model. This indicates that the linearly perturbed model (see (8)) is useful only for small time intervals.

Now let us linearize the system around the equilibrium point $E^{*}$. Under the translation

$$u=C-C^{*},v=L-L^{*}\mathrm{and}w=P-P^{*},$$

system (5) becomes

$$dY\left(t\right)=\tilde{f}\left(Y\left(t\right)\right)dt+\sum _{i=1}^3{\tilde{\sigma }}_i\left(Y\left(t\right)\right)d{W}_i\left(t\right),$$

(7)

where $Y={(u,v,w)}^T$,

$$\tilde{f}:=\left(\begin{array}{c}\alpha (u+C^{*})(w+P^{*})-{\gamma }_1(u+C^{*})(v+L^{*})\\ -\beta (v+L^{*})+{\gamma }_2(u+C^{*})(v+L^{*})\\ \delta (w+P^{*})(\xi -(w+P^{*}))-ϵ(u+C^{*})\end{array}\right)$$

and

$${\tilde{\sigma }}_1=\left(\begin{array}{c}-{\sigma }_1(u+C^{*})(v+L^{*})\\ 0\\ 0\end{array}\right),{\tilde{\sigma }}_2=\left(\begin{array}{c}0\\ {\sigma }_2(u+C^{*})(v+L^{*})\\ 0\end{array}\right),{\tilde{\sigma }}_3=\left(\begin{array}{c}0\\ 0\\ -{\sigma }_3(u+C^{*})\end{array}\right).$$

Hence, the affine linear approximation of (7) is

$$dZ\left(t\right)=AZ\left(t\right)dt+\sum _{i=1}^3(b_i+B_{i}Z\left(t\right))d{W}_i\left(t\right),$$

(8)

where

$$A=\left(\begin{array}{ccc}\alpha P^{*}-{\gamma }_1{L}^{*}& -{\gamma }_1{C}^{*}& \alpha C^{*}\\ {\gamma }_2{L}^{*}& -\beta +{\gamma }_2{C}^{*}& 0\\ -ϵ& 0& \delta (\xi -2P^{*})\end{array}\right)$$

(9)

and the diffusion components are

$$b_1=\left(\begin{array}{c}-{\sigma }_1{C}^{*}{L}^{*}\\ 0\\ 0\end{array}\right),b_2=\left(\begin{array}{c}0\\ {\sigma }_2{C}^{*}{L}^{*}\\ 0\end{array}\right),b_3=\left(\begin{array}{c}0\\ 0\\ -{\sigma }_3{C}^{*}\end{array}\right),$$

and

$$B_1=\left(\begin{array}{ccc}-{\sigma }_1{L}^{*}& -{\sigma }_1{C}^{*}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),B_2=\left(\begin{array}{ccc}0& 0& 0\\ {\sigma }_2{L}^{*}& {\sigma }_2{C}^{*}& 0\\ 0& 0& 0\end{array}\right),B_3=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ -{\sigma }_3& 0& 0\end{array}\right).$$

(10)

The asymptotic behavior of the mean square of the linearization of $Y\left(t\right)$ is studied next.

Theorem 2.

Consider the symmetric matrix

$${\displaystyle M=\left(\begin{array}{ccc}{\displaystyle \frac{{\alpha }^2{\gamma }_2({\sigma }_1^2+{\sigma }_2^2){\left(P^{*}\right)}^2+\beta {\gamma }_1^2{\sigma }_2^2+{\gamma }_1^2{\gamma }_2{\sigma }_3^2}{{\gamma }_1^2{\gamma }_2}}& {\displaystyle \frac{\alpha \beta ({\sigma }_1^2+{\sigma }_2^2)P^{*}+\alpha {\gamma }_2^2{P}^{*}-2\beta {\gamma }_1^2}{2{\gamma }_1{\gamma }_2}}& {\displaystyle \frac{\alpha \beta }{{\gamma }_2}-ϵ}\\ \frac{\alpha \beta ({\sigma }_1^2+{\sigma }_2^2)P^{*}+\alpha {\gamma }_2^2{P}^{*}-2\beta {\gamma }_1^2}{2{\gamma }_1{\gamma }_2}& \frac{{\beta }^2({\sigma }_1^2+{\sigma }_2^2)}{{\gamma }_2^2}& 0\\ {\displaystyle \frac{\alpha \beta }{{\gamma }_2}-ϵ}& 0& 2\delta (\xi -2P*)\end{array}\right)}$$

(11)

where $P^{*}$ is given in (3). Moreover, let Z be the solution of equation (8), with $Z\left(0\right)=Y\left(0\right)$, then

$${\parallel Z\left(0\right)\parallel }^2{e}^{{\lambda }_{0}t}+{\displaystyle \frac{K}{{\lambda }_0}}\left(e^{{\lambda }_{0}t}-1\right)\le \mathbb{E}\left[{\parallel Z\left(t\right)\parallel }^2\right]\le {\parallel Z\left(0\right)\parallel }^2{e}^{{\lambda }_{1}t}+{\displaystyle \frac{K}{{\lambda }_1}}\left(e^{{\lambda }_{1}t}-1\right),t>0,$$

(12)

where K is a positive constant, and ${\lambda }_0$ and ${\lambda }_1$ are the minimum and maximum eigenvalues, respectively, of the matrix M.

Proof. 

Let $V\left(t\right)=\mathbb{E}\left[{\parallel Z\left(t\right)\parallel }^2\right]$, $t>0$. Applying Itô’s formula to $Z^{\top }Z$ we obtain

$$d\left(Z^{\top }Z\right)=2Z^{\top }dZ+\sum _{i=1}^3{\parallel B_{i}Z+b_i\parallel }^{2}dt$$

and by taking expectation,

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dV}{dt}}& =& \mathbb{E}\left[Z^{\top }(A^{\top }+A)Z\right]+\sum _{i=1}^3\mathbb{E}\left[\parallel B_i{Z\parallel }^2\right]+\sum _{i=1}^3{\parallel b_i\parallel }^2\hfill \\ & =& \mathbb{E}\left[Z^{\top }\left(A^{\top }+A+\sum _{i=1}^3{B}_i^{\top }{B}_i\right)Z\right]+\sum _{i=1}^3{\parallel b_i\parallel }^2.\hfill \end{array}$$

(13)

Let $M:=A^{\top }+A+{\sum }_{i=1}^3{B}_i^{\top }{B}_i$. The Rayleigh inequality implies

$${\lambda }_0\mathbb{E}\left[Z^{\top }Z\right]\le \mathbb{E}\left[Z^{\top }MZ\right]\le {\lambda }_1\mathbb{E}\left[Z^{\top }Z\right].$$

From the inequality above and equality (13) it follows that

$${\lambda }_{0}V+K\le {\displaystyle \frac{dV}{dt}}\le {\lambda }_{1}V+K,$$

where $K:={\sum }_{i=1}^3{\parallel b_i\parallel }^2$. Hence

$${\displaystyle \frac{d}{dt}}\left(\left(V+{\displaystyle \frac{K}{{\lambda }_1}}\right)e^{-{\lambda }_{1}t}\right)\le 0\le {\displaystyle \frac{d}{dt}}\left(\left(V+{\displaystyle \frac{K}{{\lambda }_0}}\right)e^{-{\lambda }_{0}t}\right),$$

which implies

$$\left(V+{\displaystyle \frac{K}{{\lambda }_1}}\right)e^{-{\lambda }_{1}t}-\left(V\left(0\right)+{\displaystyle \frac{K}{{\lambda }_1}}\right)\le 0\le \left(V+{\displaystyle \frac{K}{{\lambda }_0}}\right)e^{-{\lambda }_{0}t}-\left(V\left(0\right)+{\displaystyle \frac{K}{{\lambda }_0}}\right),$$

from which inequality (12) follows easily.

On the other hand, using the definitions (9) and (10) of the matrices A and $B_i$, respectively, we obtain

$$M=\left(\begin{array}{ccc}{\left(L^{*}\right)}^2({\sigma }_1^2+{\sigma }_2^2)+{\sigma }_3^2-2{\gamma }_1{L}^{*}+2\alpha P^{*}+{\sigma }_2^2{C}^{*}& L^{*}{C}^{*}({\sigma }_1^2+{\sigma }_2^2)-{\gamma }_1{C}^{*}+{\gamma }_2{L}^{*}& \alpha C^{*}-ϵ\\ \\ L^{*}{C}^{*}({\sigma }_1^2+{\sigma }_2^2)-{\gamma }_1{C}^{*}+{\gamma }_2{L}^{*}& C^{*2}({\sigma }_1^2+{\sigma }_2^2)+2{\gamma }_2{C}^{*}-2\beta & 0\\ \\ \alpha C^{*}-ϵ& 0& 2\delta (\xi -2P^{*})\end{array}\right).$$

From (1) and (3) we recover version (11) of matrix M.    □

5. The Exit Time

In this section, we introduce and discuss several fundamental concepts and results. While these are well established in the literature (see, for example, [10]), our primary aim is to lay the theoretical groundwork essential for the subsequent sections.

Let $D\subset {\mathbb{R}}^3$ be a bounded domain. Consider the differential operator

$$Lu:={\displaystyle \frac{1}{2}}\sum _{i=1}^3\sum _{j=1}^3{a}_{i,j}{\displaystyle \frac{{\partial }^{2}u}{\partial x_i\partial x_j}}+\sum _{i=1}^3{f}_i{\displaystyle \frac{\partial u}{\partial x_i}},u\in C^2\left(D\right),$$

(14)

where the vector $f={(f_1,f_2,f_3)}^{\top }$ is defined in (6), and the matrix $A={\sigma }^{\top }\sigma$ is given by

$$\left(a_{ij}\right):=\left(\begin{array}{ccc}\hfill {\sigma }_1^2{C}^2{L}^2& 0& 0\hfill \\ \hfill 0& {\sigma }_2^2{C}^2{L}^2& 0\hfill \\ \hfill 0& 0& {\sigma }_3^2{C}^2\hfill \end{array}\right).$$

(15)

For the remainder of this work, our domain of interest D will be a rectangular parallelepiped, specifically a set of the form $D=(a_1,a_2)\times (b_1,b_2)\times (c_1,c_2)$. It is important to note, however, that the subsequent results apply in more general contexts. We assume that within the rectangular parallelepiped D the following condition holds:

$$\mu :=\underset{(C,L,P)\in \overline{D}}{inf}\{{\sigma }_1^2{C}^2{L}^2,{\sigma }_2^2{C}^2{L}^2,{\sigma }_3^2{C}^2\}>0,$$

where $\overline{D}$ denotes the closure of D. This condition ensures that the differential operator L is uniformly elliptic in D, implying that

$$\sum _{i=1}^3\sum _{j=1}^3{a}_{i,j}{\xi }_i{\xi }_j={\sigma }_1^2{C}^2{L}^2{\xi }_1^2+{\sigma }_2^2{C}^2{L}^2{\xi }_2^2+{\sigma }_3^2{C}^2{\xi }_3^2\ge \mu {\left|\xi \right|}^2,$$

(16)

for all $\xi ={({\xi }_1,{\xi }_2,{\xi }_3)}^{\top }\in {\mathbb{R}}^3$. Consequently (see [9]), there exists a unique function u of class $C\left(\overline{D}\right)\cap C^2\left(D\right)$ satisfying:

$$\left\{\begin{array}{cc}Lu\left(x\right)=-1,\hfill & x\in D,\\ u\left(x\right)=0,\hfill & x\in \partial D,\end{array}\right.$$

(17)

with $\partial D$ being the boundary of domain D.

Next, we introduce the mathematical object of primary interest: the first exit time of the process X from the bounded domain D, defined as

$${\tau }_D=inf\{t\ge 0:X\left(t\right)\notin D\}.$$

It is a well-established result (see [10]) that if ${inf}_{(C,L,P)\in \overline{D}}{a}_{i,i}(C,L,P)>0$, then

$${\mathbb{E}}^x\left[{\tau }_D\right]<\infty ,\mathrm{for}\mathrm{all}x\in D,$$

(18)

where the expectation ${\mathbb{E}}^x$ is taken with respect to the probability measure ${\mathbb{P}}^x[·]:=\mathbb{P}[·|X\left(0\right)=x]$.

Theorem 3.

Let u be the solution to (17). Then,

$$u\left(x\right)={\mathbb{E}}^x\left[{\tau }_D\right],\mathrm{for}\mathrm{all}x\in \overline{D}.$$

(19)

Proof. 

From (18), we know that ${\tau }_D<\infty$ almost surely. By applying Itô’s formula, the result follows. For further details, see [8].    □

For the subsequent result, it is useful to introduce the following parabolic differential operator

$$\overline{L}v:=Lv-{\displaystyle \frac{\partial v}{\partial t}},v\in C^2(D\times (0,\infty )).$$

In [7], it is proven that a unique function v of class $C(\overline{D}\times (0,\infty ))\cap C^2(D\times (0,\infty ))$ exists and satisfies:

$$\left\{\begin{array}{cc}\overline{L}v(x,t)=0,\hfill & x\in D,t>0,\\ v(x,t)=0,\hfill & x\in \partial D,t\ge 0,\\ v(x,t)=1,\hfill & x\in D,t=0.\end{array}\right.$$

(20)

With this result in hand, we can determine the distribution of the first exit time.

Theorem 4.

Let v be the solution of (20). Then,

$$v(x,t)={\mathbb{P}}^x[{\tau }_D>t],for allx\in \overline{D}andt>0.$$

(21)

Proof. 

This result also follows from Itô’s formula. For more details, see [12] or [5].    □

6. Solving the PDEs (17) and (20) Using FreeFEM

This section presents the numerical algorithm used to solve the partial differential equations (PDEs) (17) and (20), which is based on the FreeFEM software.

The algorithm implemented in FreeFEM for solving PDE (17) is structured into three main components (see [6]):

I. 

Definition of the domain D.

II. 

Definition of the functions $f_i,a_{i,i}$ and $\partial a_{i,i}/\partial x_i$.

III. 

Definition of the problem:

•    int3d(Th)$\left(G(u,w)\right)$; (Bilinear form)
•    - int3d(Th)(w); (Linear form)
•    + on(1, 2, 3, 4, 5, 6, $u=0$); (Boundary condition)

We now present the specific FreeFEM code used to solve equation (17) with the parameter values given in Section 7:

load "msh3"

// Domain definition $D=(x_1,x_2)\times (y_1,y_2)\times (z_1,z_2)$

real x1 = 0.5, x2 = 0.8, y1 = 0.15, y2 = 0.35, z1 = 0.75, z2 = 0.85;

int k = 48;    // Number of subdivisions of D

mesh3 Th = cube(k, k, k, [x*(x2-x1)+x1,y*(y2-y1)+y1,z*(z2-z1)+z1]);   // Mesh of D

// Parameter definitions

real xi = 0.7; epsi = 0.025; alfa = 0.1016; beta = 6.5525; gama1 = 0.1463;

real gama2 = 9.4913; real delta =0.15; sigma1 = 0.01118; sigma2 = 0.08131; sigma3 = 0.26811;

// Definition of coefficient functions

func a11 = sigma1*sigma1*x*x*y*y;

func a22 = sigma2*sigma2*x*x*y*y;

func a33 = sigma3*sigma3*x*x;

func da11 = 2*sigma1*sigma1*x*y*y;

func da22 = 2*sigma2*sigma2*x*x*y;

func b1 = alfa*x*z - gama1*x*y;

func b2 = - beta*y + gama2*x*y;

func b3 = delta*z*(xi-z) - epsi*x;

fespace Vh(Th, P1);    Vh u, w;    // Function space on the 3D mesh

func g = 0;    // Boundary condition

// Variational problem definition

problem Problem(u, w) = int3d(Th)(

0.5*a11*dx(u)*dx(w) + 0.5*a22*dy(u)*dy(w) + 0.5*a33*dz(u)*dz(w)

+ 0.5*da11*dx(u)*w + 0.5*da22*dy(u)*w

- b1*dx(u)*w - b2*dy(u)*w - b3*dz(u)*w)    // Bilinear part

- int3d(Th)(w)    // Linear part

+ on(1, 2, 3, 4, 5, 6, u = g);    // Dirichlet boundary conditions

Problem;    // Solve the problem

ofstream fout("CorruE.txt");    // Open output file

real x0 = 0.69; // Section for plotting

real epsilon = 1e-2;    // Tolerance for comparison

// Evaluate u on the plane $x\approx x_0$

for (int i = 0; i < Th.nv; i++)

// Save only the points where $x\approx x_0$

if (abs(Th(i).x - x0) < epsilon) {

fout $<<$ Th(i).y $<<$ " " $<<$ Th(i).z $<<$ " " $<<$ u(Th(i).x, Th(i).y, Th(i).z) $<<$ endl;

// Store $(y,z,u)$    }    }

To solve the parabolic PDE (20), the time interval is discretized, and the corresponding elliptic problem is solved at each step. Consequently, to steps $\mathbf{I}$–$\mathbf{III}$ above, a fourth step must be added to the FreeFEM scheme:

IV. 

Iterative loop. Solve the elliptic PDE, applying steps $\mathbf{I}$–$\mathbf{III}$, over the time partition.

Since the resulting code is structurally similar to the previous case, we omit the explicit FreeFEM implementation for (20). Interested readers may consult [5] and [6] for further details on the construction of such numerical schemes.

7. A Stochastic Model for Corruption in Mexico

In this section, we study the stochastic models developed earlier applied to a real-world problem, namely corruption in Mexico. Using data from Transparency International and the Global Corruption Index, we compiled Table 1, which contains information on the perception of corruption in Mexico (see also [14], where some values are estimated).

As a starting point, we use the information from the year 2022, that is, $\left(C\right(0),L(0),P(0\left)\right)=(0.69,0.3992,0.8)$. We set $P\left(0\right)=0.8$ under the assumption that current anti-corruption policies are somewhat lax (see, for example, [14]). On the other hand, we know that the parameter $\delta$ in the deterministic case controls the speed of stabilization of the laxity of anti-corruption policies. In view of this, we propose (in concordance with [14]) a relatively small value, $\delta =0.15$. The rationale is that policies and social practices generally require some time to be accepted and enforced by the population.

From the data in Table 1, we also see that the normalized corruption perception index varies very little. Similarly, there is little variability in the proportion of corrupt acts that are reported and sanctioned. This suggests taking the population standard deviations for the parameters ${\sigma }_1$ and ${\sigma }_2$, thus ${\sigma }_1=0.01118$ and ${\sigma }_2=0.08131$. Our interest lies in the behavior of the stochastic model near the unique asymptotically stable equilibrium of the deterministic system. As argued, this is reasonable for a political system, where stability is expected at least in the medium term, even though some random perturbations occur. Accordingly, we propose values for ${\sigma }_3$ considering the population standard deviation of $P\left(0\right)$ and $P^{*}$. Moreover, following [14], we consider the values $0.7$, $0.4$ for the parameter $\xi$ and $0.025$, $0.005$ for $ϵ$. The parameter values are summarized in Table 2.

Using the initial condition $Z\left(0\right)=(C\left(0\right)-C^{*},L\left(0\right)-L^{*},P\left(0\right)-P^{*})$, Theorem 2 gives

$$L\left(t\right)\le \mathbb{E}\left[{\parallel Z\left(t\right)\parallel }^2\right]\le U\left(t\right),t\ge 0,$$

where

$$L\left(t\right):=A_0{e}^{{\lambda }_{0}t}+B_0,\mathrm{and}U\left(t\right):=A_1{e}^{{\lambda }_{1}t}-B_1.$$

The parameter values in L and U are provided in Table 3 for Rows 1, 2, and 3 of Table 2. Row 4 is not considered since no equilibrium point exists for those parameters.

As is well known, linear approximations often provide only local information, but we aim to make the most of it. From Figure 1, Figure 2 and Figure 3 we obtain Table 4, Table 5 and Table 6, respectively. These show that in mean square the linear approximation has greater dispersion when the parameters of Row 2 are used in the model. Conversely, we observe that the most stable model is obtained for the parameter values in Row 1. This behavior is not at all predictable. Indeed, the parameter values ${\sigma }_1$ and ${\sigma }_2$ are significantly smaller than ${\sigma }_3$ (see Table 2); therefore, the stochastic behavior (and consequently the mean square error) would be expected to be driven mainly by this parameter. What actually occurs, however, is the opposite: the smaller ${\sigma }_3$ is, the larger the “error.”

The resolution of this apparent “paradox’’ lies in a deeper understanding of the dynamics of the linearized system. The key factor determining the exponential growth rate of the error is not the direct magnitude of the dispersion parameter ${\sigma }_3$, but rather the eigenvalues $({\lambda }_0,{\lambda }_1)$ of the matrix M, which, besides ${\sigma }_3$, depends on other parameters. Thus, an increase in ${\sigma }_3$ does not necessarily lead to a significant change in the eigenvalues of M. This result highlights a fundamental fact about complex dynamical systems: their behavior and stability are not linear functions of a single parameter. Instead, they are emergent, collective properties of the nonlinear interplay among all parameters. A small variation in one variable can significantly alter the eigenvalue structure of the system, generating effects that are not intuitively predictable, as in this case.

In what follows, we analyze row by row the parameter values from Table 2.

Row 1 

(Stability): With a relatively large $\xi$ (0.7) and a small $ϵ$ (0.005), the government’s response to corruption is weak. This parameter combination yields an equilibrium point $P^{*}=0.0345$, which reflects very strict policies. The associated matrix M has a maximum eigenvalue (${\lambda }_1$) of 0.2563, the lowest among the three scenarios. As a result, the difference $U\left(t\right)-L\left(t\right)$ grows very slowly and in a controlled manner, as shown in Table 4. In other words, the linearized model remains valid only for a very short period of time.

Row 3 

(Moderate stability): Reducing the maximum permissiveness to $\xi$$\left(0.4\right)$ while keeping $ϵ$ small (0.005) yields a similarly low equilibrium level of laxity $P^{*}$$\left(0.0696\right($. Although the volatility parameter ${\sigma }_3$ is smaller than in Row 1 (0.36520 vs 0.38275), the corresponding matrix M produces a maximum eigenvalue ${\lambda }_1$$\left(0.4704\right)$, which is larger than in Row 1. Therefore, the error growth in Table 6 is faster than in Table 4, but still much slower than in Row 2.

Row 2 

(High instability): This is the most revealing scenario. With a large $\xi$ (0.7) and a large $ϵ$ (0.025), the government’s response is much stronger than in the other cases. However, this combination of parameters leads to an equilibrium point with significantly higher laxity $P^{*}$ (0.2637). Despite having the lowest value of ${\sigma }_3$ (0.26811), this particular equilibrium configuration produces a matrix M with an exceptionally high maximum eigenvalue (${\lambda }_1$) of 1.6748. It is precisely this disproportionate eigenvalue that drives the exponential growth of the error, as shown in Table 5, where the difference $U\left(t\right)-L\left(t\right)$ grows explosively over time.

From a public policy perspective, this result suggests that merely intensifying the institutional response to perceived corruption is not sufficient: an overly aggressive response without effective sanctions may increase the observed uncertainty. Therefore, prioritizing the effectiveness of sanctions (${\gamma }_1$) and institutional robustness (${\gamma }_2$), while designing graded responses (controlling $\epsilon$), is more effective in stabilizing the dynamics than simply increasing the reaction intensity. Finally, it should be emphasized that the above conclusions stem from the linear approximation around the equilibrium $E^{*}$ and are therefore valid only locally in time; it is recommended to complement these findings with sensitivity analyses (parameter sweeps in $\epsilon ,{\gamma }_1,{\sigma }_3$) and nonlinear simulations to test the robustness of the policy recommendations outside the linear regime.

To study the global behavior of the stochastic system, we now examine the exit times from the following domain:

$$D=(0.1,0.9)\times (0.1,0.9)\times (0.1,0.9).$$

(22)

Since the function $(C,L,P)\mapsto E^{(C,L,P)}\left[{\tau }_D\right]$ cannot be easily plotted, and noting that the parameter that changes the least over time is C, we fix $C\left(0\right)=0.69$ and consider the section $D_{0.69}$ of D to obtain information about the mean exit time. Table 7 reports, for each row of parameters from Table 2, the mean exit time $u_i$ when starting at $x_0={(C\left(0\right),L\left(0\right),P\left(0\right))}^{\top }={(0.69,0.3992,0.8)}^{\top }$. Moreover, $u_{max}$ represents the point $(L,P,u)$ where the expected exit time $u=E^{(0.69,L,P)}\left[{\tau }_D\right]$ attains its maximum in $D_{0.69}$.

From Table 7 we observe that the mean exit time is greatest in Row 2, which exhibits the lowest volatility; in this case, the average exit time is $2.46$ time units. Moreover, the maximum expected exit time is obtained at the point $(0.69,0.11,0.53)$. This is noteworthy, as it indicates that when the legal system L tends to tolerate impunity and anticorruption policies are relatively lax, the system takes longer to exit its boundaries, i.e., to leave D, see (22). The longest possible residence time is $6.13$ time units.

Figure 4, Figure 5, and Figure 6 display the plots of the functions $(L,P)\mapsto E^{(0.69,L,P)}\left[{\tau }_D\right]$ and $t\mapsto P^{x_0}[{\tau }_D>t]$, respectively. Using this information, we fill Table 8 with the values of $P^{x_0}[{\tau }_D>t]$ for different times and for the three parameter rows of Table 2.

From Table 8 we observe that, for short times, the most stable rate of change of $P^{x_0}[{\tau }_D\ge t]$ corresponds to Row 2, whereas the greatest variability occurs in Row 1, particularly for $t\le 0.5$. On the other hand, for longer times ($t\approx 0.6$), there is no monotonic relation, as in the previous cases, with respect to ${\sigma }_3$. Specifically, with lower uncertainty, it is more likely that the dynamics have already left the domain D, i.e., the process exits faster ($P^{x_0}[{\tau }_D\ge 0.6]\approx 0.52$). Conversely, under higher uncertainty—as in Row 1 ($P^{x_0}[{\tau }_D\ge 0.6]\approx 0.64$)—the probability of remaining within the domain for a longer time increases.

The analysis carried out above depends only on the initial point $x_0$. However, from the first part of Figure 4, Figure 5 and Figure 6, we observe the global behavior of $(L,P)\mapsto u=E^{(0.69,L,P)}\left[{\tau }_D\right]$. Here we see that the data from Row 1 exhibit a more dramatic behavior: near the boundaries the function u is almost 0, and then it increases rapidly. Moreover, we observe that along the P-sections the values of u remain essentially constant in L, meaning that the influence of L is relatively marginal. This implies that the expected average exit time from D is not significantly affected when the perception of corruption and the laxity of anticorruption policies are fixed.

The evidence also shows that the variable that changes the least over time is the normalized index of corruption perception. Next comes the proportion of corrupt acts that have been prosecuted and sanctioned by law, and finally, the laxity of anticorruption policies. To confirm this numerically, let us consider the following sets:

$$\begin{array}{c}{D}_1=(0.6,0.8)\times (0.3,0.5)\times (0.7,0.9),\\ D_2=(0.6,0.8)\times \left(\mathbf{0}.\mathbf{2},\mathbf{0}.\mathbf{6}\right)\times (0.7,0.9),\\ D_3=(0.6,0.8)\times (0.3,0.5)\times \left(\mathbf{0}.\mathbf{6},\mathbf{1}\right).\end{array}$$

In sets $D_2$ and $D_3$ we enlarge by $0.1$ the margin of set $D_1$ in the variables L and P, respectively. Recall that the initial point is $x_0={(0.69,0.3992,0.8)}^{\top }$. Table 9 shows several expected values of exit times $u_i=E^{x_0}\left[{\tau }_{D_i}\right]$ for the parameters of Row i.

From Table 9 we confirm that random perturbations are more pronounced in the variable P: even a small change in the boundaries, from $(0.7,0.9)$ to $(0.6,1)$, produces significant variations in the mean exit time, and this holds for all three parameter rows. Furthermore, the same table shows that this behavior is indeed related to the value of ${\sigma }_3$: the smaller the value of ${\sigma }_3$, the larger the expected exit time. For this model, this indicates that greater volatility leads to a faster exit from the domain of interest, which in this case agrees with intuition. In terms of the corruption model, we may conclude that the laxity of anticorruption policies is strongly affected if the normalized index of corruption perception is modified.

It is worth noting that there is no contradiction with the observations from the linearized model, since both reflect different aspects of the dynamics. The linear approximation around $E^{*}$ evaluates local stability (growth of the mean squared error) and is dominated by the eigenvalues of matrix M, which depend on the interaction of all parameters. In contrast, the exit-time analysis from D examines the global behavior under noise: greater volatility in ${\sigma }_3$ induces stronger fluctuations that, on average, drive the system out of the domain more quickly. Hence, a larger ${\sigma }_3$ accelerates the global exit without necessarily implying a proportional increase in the eigenvalues of M (and thus in local instability). In summary, local stability around $E^{*}$ (governed by the eigenvalues of M) and the global exit speed from D (influenced by ${\sigma }_3$) are distinct but coherent aspects of the model.

From a public policy perspective, these findings highlight the importance of distinguishing between local stability and global resilience of the system. A context of high volatility in the laxity of anticorruption policies (${\sigma }_3$ large) implies that, although the system may appear locally stable around the equilibrium $E^{*}$, in practice it will exit the admissible domain D more quickly, leading to recurrent fluctuations in governance quality. Thus, effective policy design should not only aim at reducing immediate instability (through adjustments in $\epsilon$ or ${\gamma }_1$) but also at containing long-term volatility by strengthening institutional capacity to sustain consistent enforcement. In other words, anticorruption strategies should combine robust sanction mechanisms with mechanisms that minimize the dispersion in policy responses, thereby ensuring both local stability and global durability of the system’s dynamics.

8. Conclusions

The linearized analysis of the stochastic model around the deterministic equilibrium $E^{*}$ shows that the growth of the system’s second moment (variance) is critically determined by the eigenvalues of the linearized Jacobian matrix. In particular, the presence of a dominant eigenvalue with a positive real part leads to exponential growth of variability, thereby limiting the temporal validity of the linear approximation (quadratic bounds). This highlights that even small changes in institutional or punitive parameters can unpredictably amplify the dispersion of the system’s dynamics.

The numerical application to corruption indices in Mexico further illustrates that the exit times ${\tau }_D$ (the first instance at which corruption falls below a critical threshold) depend nonlinearly on the interplay between institutional laxity (P), the effectiveness of punishment (${\gamma }_1,{\gamma }_2$), and governmental responsiveness ($ϵ$). For example, scenarios characterized by high laxity (P large) and low effective sanctions (L small) exhibit the longest mean exit times ($\approx 6.13$ units). Institutional volatility (${\sigma }_3$) also plays a complex role: while it may accelerate recovery over short horizons, its long-term effects are non-monotonic—very aggressive interventions without effective sanctions can increase uncertainty. From a public policy standpoint, these results suggest that interventions should account for the full spectral structure of the system rather than focusing on a single parameter, combining improvements in punishment effectiveness (${\gamma }_1,{\gamma }_2$) with gradual institutional responses (controlling $ϵ$) to avoid overcorrections.