Let me prefix this by saying that, yes, Cox-Ingersoll-Ross (C.I.R.) is deprecated when used to model interest rates. Yet integrals of the form $$P(0,T) = E\left(\exp\left(-\int_0^Tr_s ds\right)\right) = E\left(\exp\left(-\int_0^Tr_s ds\right)| r_0=y\right)$$ are still of interest in fields like credit risk. See for instance the paper "Affine term structure models: A time-change approach with perfect fit to market curves" by Mbaye and Vrins from 2022 .
Now if the dynamics of $r$ are given by $$dr_t= \kappa(\beta-r_t)dt + \sigma\sqrt{r_t} dW_t.$$ then one can derive $$P(0,T) = \exp\left(-\int_0^T f(0,t) dt\right)$$ with instantaneous forward curve $$f(0,t) = \frac{2 \kappa\beta\left(e^{t\gamma}-1\right)}{2\gamma+(\kappa+\gamma)\left(e^{t\gamma}-1\right)}+r_0\frac{4\gamma^2}{\left(2\gamma+(\kappa+\gamma)\left(e^{t\gamma}-1\right)\right)^2}$$ in which $\gamma$ denotes the constant $\sqrt{2\sigma^2+\kappa^2}$. This can be seen in section 4.2.2 of the linked paper or by derivation using Feynman-Kac. Now we get the following asymptotic behaviour for $t\rightarrow \infty$ given by $$\lim_{t\rightarrow \infty}f(0,t)=\frac{2\kappa\beta}{\kappa+\gamma}.$$ At the same time $r$ converges in law to a Gamma distribution $\Gamma(\alpha, \beta)$ with $\alpha=\frac{2\kappa \beta}{\sigma^2}$ and $\beta=\frac{2\kappa}{\sigma^2}$. In particular this implies, that $r_\infty$ is $\beta$ in the mean. Even worse, when simulating paths of $r$ and computing $P(0,T)$ using Monte Carlo $$\overline{P}(0,T) = \sum_{n=1}^{N_1}\frac{\exp\left(-\sum_{k=1}^{N_2}\frac{r_{\frac{kT}{N_2}}+r_{\frac{(k-1)T}{N_2}}}{2}\right)}{N_1}$$ the resulting value for some large $T$ is off in the exponent by exactly $\frac{\kappa+\gamma}{2 \kappa}$, i.e. $\overline{P}(0,T)\sim\left(P(0,T)\right)^{\frac{\kappa+\gamma}{2 \kappa}}$.
For illustration purposes, let's consider a numerical example,$r_0=0.001$, $\beta = 1$, $\sigma=5$ and $\kappa = 0.8$, $N_1=10000$ and $N_2 = 1000$, then the first graph, compares the analytic forward curve (in blue) to the one obtained by a Euler-Maruyama (in orange). I use Euler-Maruyama, despite the state transition be known to be chi-squared distributed, because I'm actually interested in a extended C.I.R., the distribution of which is unknown. First note, that the high volatility is barely visible, due to the large choice of $N_1$. Despite this, the two curves differ in both the asymptotic mean and the speed of convergence to that mean.
Finally the last graph compares the analytically obtained bond prices to the values $\overline{P}(0,T)$. Again they are miles apart.
Now it's not that surprising, that we pay a premium for higher volatility. Afterall due to convexity of the exponential function Jensen's inequality implies $E(e^{X})\geq e^{E(X)}$ and in the case of Gaussian $X$, we have $E(e^{X}) = e^{E(X)+\frac{Var(X)}{2}} > e^{E(X)}$ (see for instance Monte Carlo, convexity and Risk-Neutral ZCB Pricing). Yet in the case of Gaussian $X$ I can still use direct simulation in the way described above. And this is the question here: How can I use Monte Carlo to price Zero Coupon Bonds in the Cox-Ingersoll-Ross model, if the naive way does not work?
One more thing to note here: First I can use a change of measure $\frac{dQ_T}{dP}=\exp\left(-\frac{\gamma-\kappa}{2\sigma}\int_0^T\sqrt{r_t} dt\right)$. Using Girsanov the SDE can be transformed to $$dr_t = \frac{\kappa+\gamma}{2}\left(\frac{2\kappa \theta}{\kappa+\gamma}-r_t\right)dt +\sigma\sqrt{r_t}dW^{Q}_t. $$ If I repeat the Monte Carlo with $\kappa' = \frac{\kappa+\gamma}{2}$ and $\theta' = \frac{2\kappa\theta}{\kappa+\gamma}$, I get the following forward rate and bond prices
While the change of measure, not only adjusts for the asymptotic mean, but also the speed of mean reversion, it can be observed, that the rate of mean reversion is slightly higher in the deterministic forward curve. Furthermore even if applicable, it's completely unclear, how this measure $Q$ is motivated.
