Why a Simulated Hull–White Path Cannot Predict a Bond Yield

# Using the Hull-White model to find the 2026 value in the US bond plot


# Using the Hull-White model to find the 2026 value in the US bond plot












> Background Specifically, the 2026 value of $\sim5.01\%$

### Deriving the model

We can derive the Hull–White analytical solution

$$r_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:$$dr_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:$$dr_t + \kappa r_t dt = \kappa \theta dt + \sigma dW_t$$Multiply both sides by the integrating factor $e^{\kappa t}$ to simplify the left-hand side:

$$e^{\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^{\kappa t} r_t$:$$d\left( e^{\kappa t} r_t \right) = \kappa \theta e^{\kappa t} dt + \sigma e^{\kappa t} dW_t$$

Step 2 (Integrate from $s$ to $t$)

Integrate both sides from a starting time $s$ to a future time $t$:$$\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^{\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 $r_t$)

Divide the entire equation by $e^{\kappa t}$ (which is equivalent to multiplying by $e^{-\kappa(t-s)}$):$$r_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):$$\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)$ grows large, $e^{-\kappa(t-s)} \to 0$, meaning the expected future rate converges entirely to the long-term mean $\theta$.

- If current rates are high ($r_s > \theta$), the exponential decay pulls them down. If they are low ($r_s < \theta$), they are pulled up.

Stochastic Component (Random Shocks):$$\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:$$\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 = 18$ years, starting from $r_0 = 4.65\%$ in 2008).

Our parameters:

- Initial yield ($r_s$): $4.65\%$

- Speed of mean reversion ($\kappa$): $0.35$

- Volatility ($\sigma$): $0.55$

- Time horizon ($t - s$): $18$ years

- For the long-term mean, because it shifted from $3.0$ to $4.2$ post-2021, the effective long-term target over the latter half averages out roughly near $\theta \approx 4.2\%$.

Step 1 (Calculate the Exponential Decay Factor)$$e^{-\kappa(t-s)} = e^{-0.35 \times 18} = e^{-6.3} \approx 0.00183$$(Notice that because $18$ years is very long relative to $\kappa = 0.35$, this term becomes almost zero, meaning the starting rate $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 ($\sigma_{\text{total}}$):$$\sigma_{\text{total}} = \sqrt{0.4321} \approx 0.657\%$$

Step 4 (Final Simulated Outcome)

In the python code below, the random shock drawn ($Z \sim \mathcal{N}(0,1)$) happened to be positive (about $+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: ​



- Expected Deterministic Mean: $4.20\%$ (the baseline long-term target under the Hull–White regime-shift model)

- Total Standard Deviation ($σ$ total): $0.657\%$ (the cumulative volatility calculated via Itô's isometry over the $18$-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 correct?

I believe that I have an error with SDE used $r = 4.65, \theta = 4.2$, and $\sigma = 0.55$, apparently in percentage points. That's internally possible, but then the reported standard deviation calculation is wrong in its interpretation:

$$\sqrt{\frac{(0.55)^2}{2(0.35)}}\approx 0.657$$

meaning $0.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

$$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

$$dr_t=[\theta(t)−ar_t]dt+\sigma dW_t,$$

with a time-dependent drift chosen to fit today's initial yield curve?

## Answer by almost_surely_ (score 2, accepted)

https://quant.stackexchange.com/a/85805

### Short answer

The algebra is correct — but it's the Vasicek solution, not Hull–White, and you already spotted that. Both of your self-diagnoses are right. The problem is that there are two larger issues underneath them, and the second one is the reason the answer "works".

The headline: 5.01% is not a solution, it's a seed. And a 30-year yield is not a short rate.

### 1. Your two diagnoses are both correct

Units. Yes. With $r$ in percentage points, $\sigma = 0.55$ means an annualised short-rate volatility of 55 bp, and $\sigma/\sqrt{2\kappa} = 0.657$ is 65.7 bp. That is internally consistent. Note the Feller check in your CIR block is also fine and, reassuringly, scale-invariant: $2\kappa\theta = 2.64 \ge \sigma^2 = 0.1225$ in percentage points, and $0.0264 \ge 0.001225$ in decimals. No problem there.

Vasicek vs Hull–White. Also correct, and more consequential than it looks. Hull–White (extended Vasicek) is

$$dr_t = [\theta(t) - a\,r_t]\,dt + \sigma\,dW_t,$$

where $\theta(t)$ is not a free parameter — it is pinned by today's initial forward curve $f(0,t)$:

$$\theta(t) = \frac{\partial f(0,t)}{\partial t} + a\,f(0,t) + \frac{\sigma^2}{2a}\left(1 - e^{-2at}\right).$$

That is the entire point of the model: it reprices the initial term structure exactly, by construction. Your piecewise-constant $\mu$ (3.0 before 2021, 4.2 after) is not a curve fit — it's a level chosen in hindsight because we know rates rose after 2021. So the model isn't calibrated to a curve; it's calibrated to the answer.

One consequence worth flagging: in Hull–White, $\theta(t)$ is fixed at $t=0$ from the curve you observe then. Re-picking it mid-path because of what subsequently happened isn't a Hull–White extension, it's hindsight.

### 2. 5.01% is the 94th percentile of a random number generator

This is the main thing. Your terminal value is one realisation under `np.random.seed(42)`. I reran your code across 2,000 seeds, changing nothing else:

|  | Terminal 2026 value |
| mean | 3.98% |
| sd | 0.66% |
| 5th pct | 2.91% |
| median | 3.99% |
| 95th pct | 5.06% |
| min / max | 1.17% / 6.58% |

Seed 42 lands at the 94.2nd percentile. Only 5.9% of seeds produce anything above 5.0%. Pick seed 0 or seed 7 and the "2026 solution" is somewhere around 4%, and the chart annotation would have read differently.

So the sentence to be careful with is "the 2026 solution is ~5.01%". The model doesn't produce a solution; it produces a distribution, whose centre is about 3.99% with a 90% interval of roughly [2.9%, 5.1%]. The realised 5.01% sits near the top of that. That is a legitimate and interesting statement. "The model finds 5.01%" is not.

(The plot's own source line reads "Hull–White Stochastic Simulation Model (Calibrated to Historical Regimes)", which suggests the orange path is the simulation output. If so, reading the 2026 value off it and then deriving that value from the model is circular — you'd be recovering the input.)

### 3. The $Z \approx +1.23$ story is backwards, and the mean is wrong

Two separate problems.

The mean is wrong. You computed $\mathbb{E}[r_t]$ using $\theta = 4.2$ for the whole 18 years. But your code uses $\mu = 3.0$ for the first 13 years and only then switches. Applying the mean recursion piecewise:

$$m_{2021} = 4.65\,e^{-0.35 \times 13} + 3.0\left(1 - e^{-0.35 \times 13}\right) = 3.017\%$$ $$m_{2026} = 3.017\,e^{-0.35 \times 5} + 4.2\left(1 - e^{-0.35 \times 5}\right) = \mathbf{3.994\%}$$

which matches the 2,000-seed simulated mean of 3.98% (the small gap is Euler bias — the exact transition scheme gives 3.996%, so your discretisation is fine). Your 4.20% is the mean of a different model, one where the high regime applied from 2008. So the correct $Z$ is $(5.01 - 3.99)/0.657 \approx \mathbf{1.54}$, not 1.23.

The narrative is circular. You describe $Z$ as positive "due to the aggressive 2022–2026 inflation/rate hike cycle". But `np.random.randn()` has no knowledge of the inflation cycle. $Z$ is whatever the pseudo-random stream summed to under seed 42. Computing $Z = (\text{output} - \text{mean})/\text{sd}$ and then presenting it as an explanation of the output is a restatement, not a derivation — you can always solve for the $Z$ that reproduces any number you already have. The inflation cycle enters your model only through the hand-set jump in $\mu$, and that jump is already inside the 3.99% mean.

### 4. The one that actually breaks it: a 30-year yield is not a short rate

You're modelling $r_t$ and reading it as the 30-year yield. In Hull–White those are different objects, linked by the affine bond formula:

$$P(t,T) = A(t,T)\,e^{-B(t,T)r_t}, \qquad B(t,T) = \frac{1 - e^{-a(T-t)}}{a}.$$

The continuously-compounded yield is $y(t,T) = \dfrac{-\ln A(t,T) + B(t,T)\,r_t}{T-t}$, so the sensitivity of the 30-year yield to the short rate is $B(t,T)/(T-t)$. With $a = 0.35$:

$$\frac{B(t,t+30)}{30} = \frac{(1 - e^{-10.5})/0.35}{30} = \frac{2.857}{30} \approx 0.095.$$

A 100 bp move in the short rate moves the 30-year yield by about 9.5 bp. Combined with your $\sigma$, the model implies a stationary 30-year yield volatility of roughly $0.095 \times 65.7 \approx \mathbf{6.3\ bp}$.

Your own chart shows the 30-year yield travelling from about 2.1% to about 5.3%, with moves well over 100 bp inside single years. A model implying 6 bp of annual yield vol is off by more than an order of magnitude. It only appears to fit because you never applied the bond formula — you treated $r_t$ directly as the yield, which silently sets the sensitivity to 1.0 instead of 0.095.

This is the general trap with short-rate models: fast mean reversion crushes the long end. $\kappa = 0.35$ is a half-life of $\ln 2/0.35 \approx 2.0$ years, which is why $e^{-0.35 \times 18} = 0.0018$ wipes out your 2008 starting value entirely. For the long end to move meaningfully with the short rate you need $a$ an order of magnitude smaller:

| $a$ | $B(t,t+30)/30$ | implied 30y yield sd |
| 0.35 | 0.095 | 6 bp |
| 0.10 | 0.317 | 39 bp |
| 0.05 | 0.518 | 90 bp |
| 0.01 | 0.864 | 336 bp |

In practice, desks calibrating Hull–White for long-dated exposure typically end up with $a$ in the low single-digit percent — sometimes pinned near zero or even negative — precisely because a one-factor model with meaningful mean reversion cannot simultaneously match short-rate dynamics and long-end volatility. That tension is a known limitation of the one-factor specification, not something you did wrong.

### 5. What a defensible version looks like

Pick one, and be explicit about which:

- Actually use Hull–White. Take an observed initial curve, back out $\theta(t)$ from the formula above, simulate $r_t$, and then price the 30-year bond via $P(t,T)$ to get a yield. Calibrate $a$ and $\sigma$ to swaption or cap volatilities, not by eye. You'll find $a = 0.35$ doesn't survive that.

- Model the 30-year yield directly as its own mean-reverting process. Statistically reasonable, and honest — but then it is a reduced-form yield model, not Hull–White, and it won't be arbitrage-free across maturities.

- Report a distribution, not a point. "Under this calibration the 2026 30-year yield has a median of 3.99% with a 90% interval of [2.9%, 5.1%]; the realised 5.3% sits in the upper tail" is a defensible sentence. A single simulated path labelled "solution" isn't.

The sanity check that catches all of this: rerun with 500 different seeds and see whether your conclusion survives. If the headline number moves by 200 bp when you change the seed, the headline number was never a model output.

### Verification code

```
import numpy as np

def terminal(seed):
    np.random.seed(seed)
    N, dt = 216, 1/12
    r = np.zeros(N); r[0] = 4.65
    for i in range(1, N):
        mu_t = 3.0 if 2008 + i*dt < 2021 else 4.2
        r[i] = r[i-1] + 0.35*(mu_t - r[i-1])*dt + 0.55*np.sqrt(dt)*np.random.randn()
    return r[-1]

t = np.array([terminal(s) for s in range(2000)])
print(f"seed 42        : {terminal(42):.2f}%")
print(f"mean / sd      : {t.mean():.2f}% / {t.std(ddof=1):.2f}")
print(f"5th / 95th pct : {np.percentile(t,5):.2f}% / {np.percentile(t,95):.2f}%")
print(f"percentile(42) : {(t < terminal(42)).mean()*100:.1f}")

# analytic mean under the piecewise drift
k = 0.35
m21 = 4.65*np.exp(-k*13) + 3.0*(1-np.exp(-k*13))
m26 = m21*np.exp(-k*5)  + 4.2*(1-np.exp(-k*5))
print(f"analytic mean  : {m26:.2f}%   (not 4.20%)")

# 30y yield sensitivity to the short rate
for a in (0.35, 0.10, 0.05):
    B = (1-np.exp(-a*30))/a
    print(f"a={a}: dy30/dr = {B/30:.3f}, implied 30y sd = {B/30*0.55/np.sqrt(2*a)*100:.0f} bp")
```

Output (note the CIR block in your original code consumes 215 draws before the Hull–White loop, so keep the ordering if you want to reproduce 5.01% exactly):

```
seed 42        : 4.42%
mean / sd      : 3.98% / 0.66
5th / 95th pct : 2.91% / 5.06%
percentile(42) : 74.9
analytic mean  : 3.99%   (not 4.20%)
a=0.35: dy30/dr = 0.095, implied 30y sd = 6 bp
a=0.1: dy30/dr = 0.317, implied 30y sd = 39 bp
a=0.05: dy30/dr = 0.518, implied 30y sd = 90 bp
```

Which incidentally makes the point twice over: drop the CIR block and seed 42 gives 4.42%, not 5.01%. The "solution" depends on how many random numbers you drew earlier in the script.