Background Specifically, the 2026 value of ∼5.01%\sim5.01\% Deriving the model We can derive the Hull–White analytical solution rt=rse−κ(t−s)+∫ste−κ(t−u)θ(u)du+σ∫ste−κ(t−u)dWur_t = r_s e^{-\kappa(t-s)} + \int_s^t e^{-\kappa(t-u)} \theta(u) du + \sigma \int_s^t e^{-\kappa(t-u)} dW_u where the continuous-time Vasicek SDE is given by: drt=κ(θ−rt)dt+σdWtdr_t = \kappa (\theta - r_t) dt + \sigma dW_t (Note: We use κ\kappa for the speed of mean reversion and θ\theta for the long-term mean here to follow standard mathematical notation). Step 1 (Use an Integrating Factor) Rewrite the SDE by expanding the drift term: drt+κrtdt=κθdt+σdWtdr_t + \kappa r_t dt = \kappa \theta dt + \sigma dW_t Multiply both sides by the integrating factor eκte^{\kappa t} to simplify the left-hand side: eκtdrt+κeκtrtdt=κθeκtdt+σeκtdWte^{\kappa t} dr_t + \kappa e^{\kappa t} r_t dt = \kappa \theta e^{\kappa t} dt + \sigma e^{\kappa t} dW_t Notice that the left-hand side is precisely the product rule for differentiation applied to eκtrte^{\kappa t} r_t : d(eκtrt)=κθeκtdt+σeκtdWtd\left( e^{\kappa t} r_t \right) = \kappa \theta e^{\kappa t} dt + \sigma e^{\kappa t} dW_t Step 2 (Integrate from ss to tt ) Integrate both sides from a starting time ss to a future time tt : ∫std(eκuru)=∫stκθeκudu+∫stσeκudWu\int_s^t d\left( e^{\kappa u} r_u \right) = \int_s^t \kappa \theta e^{\kappa u} du + \int_s^t \sigma e^{\kappa u} dW_u Evaluating the integrals: eκtrt−eκsrs=θ(eκt−eκs)+σ∫steκudWue^{\kappa t} r_t - e^{\kappa s} r_s = \theta \left( e^{\kappa t} - e^{\kappa s} \right) + \sigma \int_s^t e^{\kappa u} dW_u Step 3 (Solve for rtr_t ) Divide the entire equation by eκte^{\kappa t} (which is equivalent to multiplying by e−κ(t−s)e^{-\kappa(t-s)} ): rt=rse−κ(t−s)+θ(1−e−κ(t−s))+σ∫ste−κ(t−u)dWur_t = r_s e^{-\kappa(t-s)} + \theta \left( 1 - e^{-\kappa(t-s)} \right) + \sigma \int_s^t e^{-\kappa(t-u)} dW_u Deterministic Component (Expected Path): E[rt∣rs]=rse−κ(t−s)+θ(1−e−κ(t−s))\mathbb{E}[r_t \mid r_s] = r_s e^{-\kappa(t-s)} + \theta \left( 1 - e^{-\kappa(t-s)} \right) As the time horizon (t−s)(t-s) grows large, e−κ(t−s)→0e^{-\kappa(t-s)} \to 0 , meaning the expected future rate converges entirely to the long-term mean θ\theta . If current rates are high ( rs>θr_s > \theta ), the exponential decay pulls them down. If they are low ( rs<θr_s < \theta ), they are pulled up. Stochastic Component (Random Shocks): σ∫ste−κ(t−u)dWu\sigma \int_s^t e^{-\kappa(t-u)} dW_u Because it is an integral of a deterministic function with respect to a Wiener process, the stochastic term is normally distributed. Its conditional variance can be computed via Itô's isometry as: Var(rt∣rs)=σ22κ(1−e−2κ(t−s))\text{Var}(r_t \mid r_s) = \frac{\sigma^2}{2\kappa} \left( 1 - e^{-2\kappa(t-s)} \right) My attempt Let's use a 18-year time horizon ( t−s=18t - s = 18 years, starting from r0=4.65%r_0 = 4.65\% in 2008). Our parameters: Initial yield ( rsr_s ): 4.65%4.65\% Speed of mean reversion ( κ\kappa ): 0.350.35 Volatility ( σ\sigma ): 0.550.55 Time horizon ( t−st - s ): 1818 years For the long-term mean, because it shifted from 3.03.0 to 4.24.2 post-2021, the effective long-term target over the latter half averages out roughly near θ≈4.2%\theta \approx 4.2\% . Step 1 (Calculate the Exponential Decay Factor) e−κ(t−s)=e−0.35×18=e−6.3≈0.00183e^{-\kappa(t-s)} = e^{-0.35 \times 18} = e^{-6.3} \approx 0.00183 (Notice that because 1818 years is very long relative to κ=0.35\kappa = 0.35 , this term becomes almost zero, meaning the starting rate 4.65%4.65\% has virtually no impact left by 2026—the process has fully settled around the long-term mean). Step 2 (Calculate the Expected Value (Mean)) \begin{align*} \mathbb{E}[r_t \mid r_s] &= r_s e^{-\kappa(t-s)} + \theta \left( 1 - e^{-\kappa(t-s)} \right)\ &= (4.65 \times 0.00183) + 4.2 \left( 1 - 0.00183 \right)\ &\approx 0.0085 + 4.2 \left( 0.99817 \right)\ &\approx 0.0085 + 4.1923\ &\approx 4.20% \end{align*} Step 3 (Calculate the Standard Deviation) \begin{align*} \text{Var}(r_t \mid r_s) &= \frac{\sigma^2}{2\kappa} \left( 1 - e^{-2\kappa(t-s)} \right)\ &= \frac{0.55^2}{2 \times 0.35} \left( 1 - e^{-2(0.35 \times 18)} \right)\ &= \frac{0.3025}{0.7} \left( 1 - e^{-12.6} \right)\ &\approx 0.4321 \times (1 - 0.000003)\ &\approx 0.4321 \end{align*} Taking the square root gives the standard deviation ( σtotal\sigma_{\text{total}} ): σtotal=0.4321≈0.657%\sigma_{\text{total}} = \sqrt{0.4321} \approx 0.657\% Step 4 (Final Simulated Outcome) In the python code below, the random shock drawn ( Z∼N(0,1)Z \sim \mathcal{N}(0,1) ) happened to be positive (about +1.23+1.23 standard deviations above the mean due to the aggressive 2022–2026 inflation/rate hike cycle): \begin{align*} Z&=\frac{\text{Actual Final Value}\text{−Expected Deterministic Mean}}{\text{Total Standard Deviation}}\ &= \frac{5.01 - 4.20}{0.657}\ &= \frac{0.81}{0.657}\ &\approx +1.23 \end{align*} where: ​ Actual Final Value: 5.01%5.01\% (the terminal yield produced by the stochastic simulation) Expected Deterministic Mean: 4.20%4.20\% (the baseline long-term target under the Hull–White regime-shift model) Total Standard Deviation ( σσ total): 0.657%0.657\% (the cumulative volatility calculated via Itô's isometry over the 1818 -year period) \begin{align*} \text{Result} &= \text{Mean} + (Z \times \text{StdDev})\ &= 4.20% + (1.23 \times 0.657%)\ &\approx \boxed{5.01%} \end{align*} import numpy as np

Let's run a simulation for CIR (Cox-Ingersoll-Ross) and Hull-White models

np.newaxis np.random.seed(42)

N = 216 # Monthly steps from 2008 to 2026 dt = 1 / 12

1. CIR Model simulation: dr = theta * (mu - r_t)dt + sigma * sqrt(r_t) * dW_t

theta_cir = 0.4 mu_cir = 3.3 sigma_cir = 0.35 # Must satisfy Feller condition: 2 * theta * mu >= sigma^2 to stay positive

r_cir = np.zeros(N) r_cir[0] = 4.65

for i in range(1, N): # Ensure non-negative inside sqrt r_prev = max(0.0, r_cir[i-1]) dr = theta_cir * (mu_cir - r_prev) * dt + sigma_cir * np.sqrt(r_prev) * np.sqrt(dt) * np.random.randn() r_cir[i] = r_prev + dr

2. Hull-White Model (Time-varying mean theta(t) or drift to fit term structure)

Simplest time-varying mean formulation: theta(t) matches a shifting trend

theta_hw = 0.35 sigma_hw = 0.55

r_hw = np.zeros(N) r_hw[0] = 4.65

for i in range(1, N): t_val = 2008 + i * dt # Let the long-term mean drift higher post-2021 to capture the inflation regime shift mu_t = 3.0 if t_val < 2021 else 4.2

dr = theta_hw * (mu_t - r_hw[i-1]) * dt + sigma_hw * np.sqrt(dt) * np.random.randn()
r_hw[i] = r_hw[i-1] + dr

print(f"CIR Model 2026 Terminal Value: {r_cir[-1]:.2f}%") print(f"Hull-White (Regime-Shift) 2026 Terminal Value: {r_hw[-1]:.2f}%") ##CIR Model 2026 Terminal Value: 4.71% ##Hull-White (Regime-Shift) 2026 Terminal Value: 5.01% My question Would this derivation be succinct? I believe that I have an error with SDE used r=4.65,θ=4.2r = 4.65, \theta = 4.2 , and σ=0.55\sigma = 0.55 , apparently in percentage points. That's internally possible, but then the reported standard deviation calculation is wrong in its interpretation: (0.55)22(0.35)≈0.657\sqrt{\frac{(0.55)^2}{2(0.35)}}\approx 0.657 meaning 0.6570.657 percentage points if rates are measured in percentage points? Also, the Hull–White equation that I derived isn't actually the standard Hull–White specification we're subsequently describing, as derivation starts from drt=κ(θ−rt)dt+σdWt,dr_t=\kappa(\theta−r_t)dt+\sigma dW_t, where θ\theta is constant. That's essentially the Vasicek model, but the standard one-factor Hull–White is drt=[θ(t)−art]dt+σdWt,dr_t=[\theta(t)−ar_t]dt+\sigma dW_t, with a time-dependent drift chosen to fit today's initial yield curve?