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.
39,287 characters · 9 sections · 33 citation commands
Linear estimations of dynamic fixed effects logit models only with time effects
\allowdisplaybreaks
Honore2006611 suggest that neither the dynamic fixed effects binary choice model only with time trends nor that only with time dummies is point-identified after meticulously inspecting their calculation results from a linear programming method using probit models. However, instead, as shown from now, the logit models included in them are point-identified if five consecutive time periods are provided. Disappointingly, they overlook this significant discovery without showing the calculation results for the case of five time periods and without showing the results from using logit specifications.
The matching approach by Honore2000839 is not essential for achieving point-identification for dynamic fixed effects logit models. To show this, we propose (globally) root-{\it N} consistent estimators for both models only with time dummies and only with time trends after proving point-identification of these models.
For dynamic fixed effects logit models with strictly exogenous continuous explanatory variables and/or time effects, Kitazawa2022350 is the first to obtain the valid moment conditions, leveraging the properties of hyperbolic functions.\footnote{Honore2024, Dobronyi2021apr_arXiv, and Dano2023feb_arXiv also construct the equivalent moment conditions in different ways: the functional differencing approach proposed by Bonhomme20121337, the formulation of the full likelihood, and the method based on transition probabilities, respectively.} After the transformations of these moment conditions give birth to linear asymptotically normal and root-{\it N} consistent estimators for the transformations of parameters of interest in the models only with time effects, asymptotically normal and root-{\it N} consistent estimators are obtained for the parameters of interest without nonlinear optimization.
In retrospect, no estimators are heretofore proposed for point-identifying the parameters of interest in the dynamic fixed effects logit models without imposing some additional restrictions, except in the simple models\footnote{Chamberlain19853 shows that the model with neither explanatory variables nor time effects is point-identified. Further, Honore2000839 show that the model only with binary explanatory variables is point-identified as well. They propose root-{\it N} consistent estimators for the respective models.}: Honore2000839 propose a kernel-using estimator requiring that first-differences of the explanatory variables at a specific time period are continuously distributed with support near zero, which excludes the time effects and explanatory variables such as age and spouse income based on the age-based remuneration system, etc.\footnote{These variables could be used in studies on the determinants of female labor force participation.} Honore2024 show that the usage of the moment conditions mentioned in previous paragraph point-identifies parameters of interest for the model with strictly exogenous explanatory variables if a special dataset is used. Although Dobronyi2021apr_arXiv provide specific numerical examples of point-identification for the models with time effects, where they use the moment inequality conditions they propose in addition to the moment (equality) conditions in previous paragraph, it is unknown whether the models are generically point-identified or not.\footnote{The 2025 version of Dobronyi2021apr_arXiv, where T.M. Russell joins as one of the coauthors, definitely states in footnote 34 that it is not known whether the time trend model is generically point-identified or not when four time periods (incorporating the initial observations) are provided. No statements are found for the cases of five or more time periods. To the best of author's knowledge, no literatures and documents are available concerning proofs of point-identification of the time trend and time dummies models, although the analogy stated in Remark 8 in the supplement of the 2025 version of Dano2023feb_arXiv claims that the AR(2) (autoregressive of order 2) time trend model will be point-identified when six time periods are provided: we can fetch the supplement from Dano_HP's website.}
Although it remains an open question whether the models with strictly exogenous continuous explanatory variables are point-identified or not, this paper radically sorts out the controversy as to whether the models only with time effects are point-identified or not.
The rest of the paper is as follows: section (ref) constructs systems of simultaneous equations consisting of the transformed moment conditions to obtain solutions of the transformed parameters uniquely. Then, values of the original parameters are solved from their solutions uniquely. Section (ref) presents Monte Carlo illustrations giving evidence for asymptotic normality and root-{\it N} consistency of the estimators proposed in section (ref). Finally, section (ref) clears up the misunderstanding surmised from the results obtained by Honore2006611 and concludes the paper.
Throughout the paper, the subscripts $i$ ($i=1,\ldots,N$) and $t$ ($t=1,\ldots,T$) denote the individual and time period, respectively. It is assumed that the number of individuals $N \to \infty$, while the number of time periods $T$ is fixed. We use any five consecutive time periods within $T$ periods to estimate the models.
The model only with time dummies ($TD_t$ for $2 \leq t \leq T$) that gives birth to the binary dependent variable $y_{it}$ is as follows:
where the logit probability $p(\eta_i, y_{i,t-1}, TD_t)=\exp(\eta_i+\gamma y_{i,t-1}+TD_t)/(1+\exp(\eta_i+\gamma y_{i,t-1}+TD_t))$ with $\eta_i$ and $\gamma$ being the fixed effect and state dependence parameter respectively, with which $y_{it}=1$. $v_{it}$ is the disturbance, $y_{i1}$ is the initial value of the binary dependent variable, and $v_i^{t-1}=(v_{i1},\ldots,v_{i,t-1})$ with $v_{i1}$ being empty. We postulate that in this model, the original parameters (i.e., parameters of interest to be estimated) are $\gamma$, $\varDelta TD_t$, and $\varDelta TD_{t+1}$ for time $t$, where $\varDelta TD_t = TD_t - TD_{t-1}$.
Kitazawa2022350 proposes two types of moment conditions for the above model. We have those based on the {\it g-form} and {\it h-form} respectively as follows:
and
where putting $\delta = \exp(\gamma) - 1$, $\phi_t = \exp(\varDelta TD_t)$, and $\phi_t^{-1} = 1/\phi_t$,
By putting
we write
where the transformed parameters $\alpha_a^{(A)} = \phi_t$, $\alpha_b^{(A)} = \phi_{t+1}^{-1}$, $\alpha_c^{(A)} = \phi_t \phi_{t+1}$, $\alpha_d^{(A)} = \phi_t (1+\delta)$, $\alpha_e^{(A)} = \phi_t (1+\delta)^{-1}$, $\alpha_f^{(A)} = \phi_{t+1}^{-1} (1+\delta)^{-1}$, and $\alpha_g^{(A)} = \phi_t \phi_{t+1} (1+\delta)^{-1}$. Therefore, for time $t$, we have the following eight unconditional moment conditions linear in the seven parameters above:
where $\operatorname*{\mathbf{0}}$ is a zero vector. The derivation of $\iota_{it}$, $\iota_{it} ^{*}$, $\kappa_{it}$, and $\kappa_{it}^{*}$ is provided in appendix (ref). A direct computation with maxima suggests that the moment functions composing ((ref)) are linearly independent for $(y_{i,t-3}, \ldots , y_{i,t+1}) \in \{0,1\}^5$ and almost all values of the parameters $\gamma$, $\varDelta TD_t$, and $\varDelta TD_{t+1}$.\footnote{The discussion on the linear dependence of the moment functions is also conducted in Honore2024, Kruiniger2020oct_arXiv, Dobronyi2021apr_arXiv, and Dano2023feb_arXiv, including strictly exogenous continuous explanatory variables and longer lagged dependent variables.}
We can remove any one moment condition from the moment conditions ((ref)) to produce an unique just-identified moment estimator for the set of the seven parameters $\alpha_a^{(A)}$, $\alpha_b^{(A)}$, $\alpha_c^{(A)}$, $\alpha_d^{(A)}$, $\alpha_e^{(A)}$, $\alpha_f^{(A)}$, and $\alpha_g^{(A)}$ as long as the rank condition is satisfied in constructing the estimator; hereafter we call it the {\textsl A} estimator (see details in appendices (ref) and (ref)). Then, we can obtain some unique estimators (without numerical optimization) for the set of the original parameters $\gamma$, $\varDelta TD_t$, and $\varDelta TD_{t+1}$, utilizing some relationships among the seven parameters (see details in appendix (ref)). It should be noted that all these estimators are asymptotically normal and root-{\it N} consistent.
Of even more greater interest is the fact that using the moment conditions ((ref)) without $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \kappa_{it}]=0$ and ((ref)) without $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \iota_{it}]=0$ can give birth to unique just-identified moment estimators for the sets of the following six parameters $\alpha_a^{(A)}$, $\alpha_b^{(A)}$, $\alpha_c^{(A)}$, $\alpha_d^{(A)}$, $\alpha_f^{(A)}$, and $\alpha_g^{(A)}$ and the following six parameters $\alpha_a^{(A)}$, $\alpha_b^{(A)}$, $\alpha_c^{(A)}$, $\alpha_e^{(A)}$, $\alpha_f^{(A)}$, and $\alpha_g^{(A)}$, respectively; these estimators also count among the {\it A} estimators (see details in appendices (ref) and (ref)). Using these estimators, we can obtain asymptotically normal and root-{\it N} consistent estimators for $\gamma$, $\varDelta TD_t$, and $\varDelta TD_{t+1}$, as well (see details in appendix (ref)).
Alternatively, using
where the transformed parameters $\alpha_a^{(B)} = \phi_t^{-1}$, $\alpha_b^{(B)} = \phi_{t+1}$, $\alpha_c^{(B)} = \phi_t^{-1} \phi_{t+1}^{-1}$, $\alpha_d^{(B)} = \phi_t^{-1} (1+\delta)$, $\alpha_e^{(B)} = \phi_t^{-1} (1+\delta)^{-1}$, $\alpha_f^{(B)} = \phi_{t+1}(1+\delta)^{-1}$, and $\alpha_g^{(B)} = \phi_t^{-1} \phi_{t+1}^{-1} (1+\delta)^{-1}$, we have the following eight unconditional moment conditions linear in the seven parameters above:
The derivation of $\lambda_{it}$, $\lambda_{it} ^{*}$, $\mu_{it}$, and $\mu_{it}^{*}$ is provided in appendix (ref). Given that these moment conditions ((ref)) are a variant of ((ref)), we suggest that the moment functions composing ((ref)) are linearly independent as well as those composing ((ref)). A calculation with maxima also suggests this.
For these moment conditions ((ref)), we can obtain some asymptotically normal and root-{\it N} consistent estimatators for the the original parameters $\gamma$, $\varDelta TD_t$, and $\varDelta TD_{t+1}$ in the same manner as for the moment conditions ((ref)), where three types of the {\textsl B} estimators are defined: one using the moment conditions ((ref)) with any one moment condition removed for the set of the seven parameters $\alpha_a^{(B)}$, $\alpha_b^{(B)}$, $\alpha_c^{(B)}$, $\alpha_d^{(B)}$, $\alpha_e^{(B)}$, $\alpha_f^{(B)}$, and $\alpha_g^{(B)}$, one using ((ref)) without $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \lambda_{it}]=0$ for the set of the six parameters $\alpha_a^{(B)}$, $\alpha_b^{(B)}$, $\alpha_c^{(B)}$, $\alpha_d^{(B)}$, $\alpha_f^{(B)}$, and $\alpha_g^{(B)}$, and one using ((ref)) without $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \mu_{it}]=0$ for the set of the six parameters $\alpha_a^{(B)}$, $\alpha_b^{(B)}$, $\alpha_c^{(B)}$, $\alpha_e^{(B)}$, $\alpha_f^{(B)}$, and $\alpha_g^{(B)}$ (see details in appendices (ref), (ref), and (ref)).
In addition, using both the moment conditions available in ((ref)) or ((ref)) dated $t-1$ and some of the above estimators for the transformed parameters, we can point-identify the original parameter $\varDelta TD_{t-1}$ and obtain asymptotically normal and root-{\it N} consistent two-step estimators for it (see details in appendix (ref)).
If we estimate all the first-differences of time dummies available within $T$ periods (i.e., $\varDelta TD_3,\ldots,\varDelta TD_{T}$), we should use the aforementioned estimators for any sets of five consecutive time periods.
The model only with time trends ($\varphi (t-\tau)$ for $2 \leq t \leq T$ where $\varphi$ and $\tau$ are parameters) is as follows:
where the logit probability $p(\eta_i, y_{i,t-1}, t)=\exp(\eta_i+\gamma y_{i,t-1}+\varphi (t-\tau))/(1+\exp(\eta_i+\gamma y_{i,t-1}+\varphi (t-\tau)))$ with which $y_{it}=1$ and we inherit the notation from previous subsection other than the time trends. In this model, the original parameters are $\gamma$ and $\varphi$.
With $\phi =\phi_t = \phi_{t+1} = \exp(\varphi)$ in $\hbar U_{it}^{-}$ and $\hbar \Upsilon_{it}^{-}$, we define
where the transformed parameters $\alpha_a^{(C)} = \phi$, $\alpha_b^{(C)} = \phi^{-1}$, $\alpha_c^{(C)} = \phi^2$, $\alpha_d^{(C)} = \phi^{-2}$, $\alpha_e^{(C)} = \phi (1+\delta)$, $\alpha_f^{(C)} = \phi (1+\delta)^{-1}$, $\alpha_g^{(C)} = \phi^{-1} (1+\delta)$, and $\alpha_h^{(C)} = \phi^{-1} (1+\delta)^{-1}$. Therefore, for time $t$, we have the following eight unconditional moment conditions linear in the eight parameters above:
The derivation of $\varrho_{it}$, $\varrho_{it} ^{*}$, $\varsigma_{it}$, and $\varsigma_{it}^{*}$ is provided in appendix (ref). A calculation similar to the case regarding ((ref)) suggests the linear independence among the moment functions composing ((ref)) for almost all values of the parameters $\gamma$ and $\varphi$.
The just-identified estimator for the set of the eight parameters $\alpha_a^{(C)}$, $\alpha_b^{(C)}$, $\alpha_c^{(C)}$, $\alpha_d^{(C)}$, $\alpha_e^{(C)}$, $\alpha_f^{(C)}$, $\alpha_g^{(C)}$, and $\alpha_h^{(C)}$ is constructed using all the moment conditions ((ref)) as long as the rank condition is satisfied, similar to those in previous subsection; hereafter we call it the {\textsl C} estimator (see details in appendix (ref)). Then, we can obtain some unique estimators for the original parameters $\gamma$ and $\varphi$ (see details in appendix (ref)). All these estimators are asymptotically normal and root-{\it N} consistent.
In addition, the {\textsl A} and {\textsl B} estimators proposed in previous subsection point-identify the original parameters $\gamma$ and $\varphi$ using some relationships among the transformed parameters, as well. Obviously, their estimators are asymptotically normal and root-{\it N} consistent.
In this section, the finite sample behavior of estimators using the moment conditions proposed in previous section is investigated with Monte Carlo experiments. For the two models in previous section, the experiments are designed with the number of replications $R=2500$, the number of generated time periods $T=8$, and the cross-sectional sizes $N=(10\,\text{million}, \, 50\,\text{million}, \, 100\,\text{million})$. The first three cross-sections are discarded, while the last five ones are used for the estimations. The experiments are carried out with a program written in Julia developed by Bezanson2017julia.\footnote{The computation is performed using Julia 1.12.4.}
We present the tables displaying the Monte Carlo results, including mean (Monte Carlo mean), sd (Monte Carlo standard deviation), se (Monte Carlo mean of standard errors), bias (Monte Carlo bias), and rmse (Monte Carlo root mean squared error).
In each replication, estimates of the original parameters and those from using two-step estimations are calculated according to the illustrations in appendices (ref) and (ref), respectively.
The DGP (data generating process) is as follows:
The experiments are carried out with $\gamma=1$, $TD_1 = 0.1$, $TD_2=-0.1$, $TD_3=0.3$, $TD_4=-0.3$, $TD_5=-0.1$, $TD_6=0.3$, $TD_7=0.5$, $TD_8=0.2$, and $\sigma_\eta^2=0.5$. Accordingly, true values of the original parameters to be estimated (labelled true in the tables) are $\gamma=1$, $\Delta TD_6=0.4$, $\Delta TD_7=0.2$, and $\Delta TD_8=-0.3$.
In the experiments, the moment conditions ((ref)) with $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \kappa_{it}]=0$ excluded and ((ref)) with $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \lambda_{it}]=0$ excluded are used for just-identifying two distinct sets of the six transformed parameters, respectively. That is, we use both of the {\textsl A} estimator not using $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \kappa_{it}]=0$ and {\textsl B} estimator not using $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \lambda_{it}]=0$. The former is referred to as the {\textsl A--(3+7)} estimator, while the latter the {\textsl B--(1+5)} estimator.\footnote{Computational results using maxima reveal that both of the sets of the moment functions composing the {\textsl A--(3+7)} and {\textsl B--(1+5)} estimators are linearly independent except where $(\gamma, \, \varDelta TD_{t+1}) = (0, \, 0)$. In addition, both of the sets of the moment functions composing the {\textsl A--(1+5)} and {\textsl B--(3+7)} estimators are linearly independent except where $(\gamma, \, \varDelta TD_{t+1}) = (0, \, 0)$ as well, where the former and the latter use the moment conditions ((ref)) with $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \iota_{it}]=0$ excluded and ((ref)) with $\operatorname*{\mathrm{E}}[(1\;\; y_{i,t-3})' \; \mu_{it}]=0$ excluded respectively. Accordingly, these four sets of the moment functions are also linearly independent except where $(\gamma, \, \varphi) = (0, \, 0)$ for the case of the model with time trends. Consequently, we clarify that using these four estimators, the models with time dummies and with time trends are point-identified except for the sets of true parameter values $(\gamma, \, \varDelta TD_{t+1}) = (0, \, 0)$ and $(\gamma, \, \varphi) = (0, \, 0)$ respectively, as long as the rank conditions are satisfied (see also appendix (ref)). Because of these facts, the usage of these four estimators is encouraged for root-{\it N} consistent estimations of the models with time dummies and with time trends.}
Using the results from these two estimators, two distinct sets of estimation results for the original parameters and their variances are obtained from the relationships among $\alpha_a^{(A)}$, $\alpha_b^{(A)}$, and $\alpha_d^{(A)}$ and among $\alpha_a^{(B)}$, $\alpha_b^{(B)}$, and $\alpha_d^{(B)}$, respectively (see details in appendices (ref) and (ref)). These Monte Carlo results for the original parameters are shown in Tables (ref) and (ref), respectively.\footnote{The corresponding results for the transformed parameters are shown in Tables (ref) and (ref) in appendix (ref), respectively.}
Looking at the tables, all values of sd and rmse and almost all bias sizes decrease with the larger $N$, while values of se predict well those of sd. A calculation indicates that convergence rates of the estimators are root-{\it N}. Further, histograms and Q-Q plots generated from the estimates obtained using replication datasets depict asymptotic normality of the estimators.\footnote{The figures are available from the author upon request.}
In this case, $p(\eta_i, y_{i,t-1}, TD_t)$ and $q(\eta_i, TD_1)$ in section (ref) are replaced by $p(\eta_i, y_{i,t-1}, t)=\exp(\eta_i+\gamma y_{i,t-1}+\varphi (t-\tau))/(1+\exp(\eta_i+\gamma y_{i,t-1}+\varphi (t-\tau))$ and $q(\eta_i)=\exp(\eta_i+\varphi (1-\tau))/(1+\exp(\eta_i+\varphi (1-\tau)))$, respectively.
The experiments are carried out with $\gamma=1$, $\varphi=0.3$, $\tau=1$, and $\sigma_\eta^2=0.5$; the former two are true values of the original parameters to be estimated.
In the experiments, using the results from the {\textsl C} estimator, the original parameters and their variances are estimated from the relationships among $\alpha_a^{(C)}$ and $\alpha_e^{(C)}$ (see details in appendix (ref)). The Monte Carlo results are shown in Table (ref).\footnote{The corresponding results for the transformed parameters are shown in Table (ref) in appendix (ref).}
As well as for the case of previous subsection, all values of sd and rmse and all bias sizes decrease with the larger $N$, while values of se predict well those of sd. Root-{\it N} rate convergence and asymptotic normality of the estimator are confirmed as well as in previous subsection.
Honore2006611 suggest that no point-identification will be conducted for dynamic fixed effects binary choice models only with time effects, scrutinizing their results obtained with a linear programming method under probit specifications, although they say that the unavailability of point-identification is of little practical repercussion in light of the highly small-scale identified regions which they obtain.\footnote{Since the probit and logit functions are sililarly shaped, it is reasonably conceivable that the results from using logit specifications will be almost the same as those from using probit ones. Even so, it is a matter of regret that they were on the very verge of discovery of point-identification.} However, other than for special true values of parameters of interest and atypical datasets, point-identification is available as long as five consecutive time periods are provided and logit specifications are used, as seen from the results in previous sections.
As for the model with time dummies, in Figures 5 and 6 (where four consecutive time periods are used) in Honore2006611, horizontal sizes of the upper three identified regions are vanishingly small. In addition, the identified regions for $(\gamma, \delta_3)$ and $(\gamma, \delta_4)$ are congruent to that for $(\gamma, \delta_2)$; the former two are vertically lowered from the latter one by the sizes of $\delta_2 - \delta_3$ and $\delta_2 - \delta_4$ (i.e., $0.1$ and $0.2$), respectively. Further, the lower identified regions for $(\delta_3, \delta_2)$, $(\delta_4, \delta_2)$ and $(\delta_4, \delta_3)$ are almost the $45$-degree line segments with the intercepts being the values of $\delta_3 - \delta_2$, $\delta_4 - \delta_2$ and $\delta_4 - \delta_3$ (i.e., $0.1$, $0.2$ and $0.1$), respectively. Therefore, identified regions for $(\gamma, \delta_3 - \delta_2, \delta_4 - \delta_3)$ are almost singleton.\footnote{Note that $\gamma$, $\delta_2$, $\delta_3$, and $\delta_4$ in Honore2006611 correspond to $\gamma$, $TD_{t-1}$, $TD_t$, and $TD_{t+1}$ in section (ref), respectively. Accordingly, $\delta_3 - \delta_2$ and $\delta_4 - \delta_3$ correspond to $\varDelta TD_t$ and $\varDelta TD_{t+1}$, respectively.} Further, Figures 4 (where three consecutive time periods are used) and 5 for $\gamma=1$ are similar in shape, although sizable identified regions in the former cannot in any way be regarded as singleton.
As for the model with time trends, all nine identified regions for $(\gamma, \beta)$ in Figure 3 (where four consecutive time periods are used) in Honore2006611 are vanishingly small and almost singleton, although the identified regions in Figure 2 (where three consecutive time periods are used) are of non-negligible size.\footnote{Note that $\gamma$ and $\beta$ in Honore2006611 correspond to $\gamma$ and $\varphi$ in section (ref), respectively.}
It might be inferred from the foregoing that the identified regions obtained from the linear programming method will accurately reduce to singleton for both dynamic fixed effects logit models only with time dummies and only with time trends if five consecutive time periods are provided. This inference dovetails with the theoretical and experimental results in this paper.
The root-{\it N} consistent estimators proposed in this paper settled the problem whether dynamic fixed effects logit models only with time effects are point-identified or not. The discovery of root-{\it N} consistent estimators was in an open-and-shut manner materialized for dynamic fixed effects logit models only with time effects. Exactly, “the darkest place is under the candlestick."
{\bf Postscript 1:} To enhance the authenticity of the estimators proposed in this paper, we can conduct the Wald test of the restrictions among the transformed parameters using their estimates (see appendix (ref)).
{\bf Postscript 2:} Regarding Hahn2001913's\ (Hahn2001913) demonstration of impossibility of root-{\it N} consistent estimation for the dynamic fixed effects logit model only with time dummies when three time periods are provided, Honore2024 claim that this model is underidentified for three time periods. This claim is described in the framework of this paper as follows: with $y_{i,t-2}=0$, the two valid moment conditions $\operatorname*{\mathrm{E}}[\iota_{it}]=0$ and $\operatorname*{\mathrm{E}}[\kappa_{it}^{*}]=0$ in $\operatorname*{\mathrm{E}}[(\iota_{it}\;\; \iota_{it}^{*}\;\; \kappa_{it}\;\; \kappa_{it}^{*})'] = \operatorname*{\mathbf{0}}$ underidentify the three parameters $\delta$, $\phi_t$ and $\phi_{t+1}$, where we take into consideration that the moment functions $\iota_{it}^{*}$ and $\kappa_{it}$ are both zero.
{
The author would like to acknowledge the use of generative AI tools, specifically OpenAI ChatGPT and Google Gemini, for technical assistance in optimizing the Julia scripts used for Monte Carlo simulations and providing conceptual clarifications on mathematical methods. The author also extends the gratitude to the Stack Exchange community for providing valuable insights into programming challenges. In addition, the AI tools were utilized for linguistic refinement of the manuscript. All AI-generated suggestions and community-sourced solutions were critically reviewed and validated by the author. Finally, the brief comments from Jiaying Gu, Bo Honor\'e, Hugo Kruiniger, Jo\ ao Santos Silva, Elie Tamer, and Martin Weidner are gratefully acknowledged. }
{