I'm trying to code the Hull White stochastic volatility model using matlab and somewhere my code seems to mess up. I've coded the SABR model as well and that's working fine. When I compare prices obtained between my two models, there seems to be a linear relationship (which is what I expect if I've done the Hull White correctly). Hull White SABR price comparison

When I produce implied volatility smiles for the Hull White, it looks something like this: enter image description here

I'm using the pricing formula for the Hull White model:

enter image description here

This is the procedure I'm following for my code:

  1. Simulate Paths of sigma using the Milstein scheme for discretising an SDE
  2. Average points at each time incriment for all the simulations and compute \bar{V}.
  3. Find the Black Scholes Price
  4. Compute Implied Volatility.

It's definitely going wrong somewhere so if anyone has the time to look through my code I'd be soo grateful! thanks in advance

function [Imp_Vol, Price] = Hull_White(S0, K, r, sigma0,gamma,alpha, T, M)

%   Monte Carlo European Call Option with Local Volatility using the Milstein Scheme. 
%   S0 - Underlying price at T=0
%   K - Strike Price
%   T - Time to Maturity in Years (please input uniform decimals ie 1.5,
%   2.9)
%   r - Interest rate
%   sigma0 - Volatility
%   gamma - Variable in Hull White model
%   alpha - Variable in Hull White model
%   M - Number of Monte Carlo Simulations 
%   Output price - corrosponding option price

N=floor(T*260); % Days to maturity

dt=1/260; % incriment for descretization

sigma_path = zeros(N,M);

for j=1:M % starting the Monte Carlo Simulations

    Z=randn(N,1);  % N random numbers to approximate dWi
    sigma=zeros(N,1); % initialising the volatility path vector

    sigma(1) = sigma0;

     for k=2:N % Calculating the volatility Path using the Milstein Scheme

            sigma(k) = sigma(k-1)*(1 + sqrt(dt)*gamma*Z(k-1) + dt*(alpha + 0.5*gamma*gamma*(Z(k-1)^2-1)));
sigma_integral = mean(sigma_path')';
sigma_integral_squared = sigma_integral.^2;
dt_vector = dt*ones(N,1);

V_bar = mean(sigma_integral_squared'*dt_vector); %numerical integration

[Call, Put] = blsprice(S0, K, r, T, V_bar, 0);
Imp_Vol = blsimpv(S0, K, r, T, Call, [],[], [], []);
Price = Call;
  • $\begingroup$ You're discretisation scheme seems weird to me. As you've defined it, the SDE for $\sigma_t$ is a simple GBM. So why would you bother with a Milstein scheme? Euler on log-returns should be fine (ie replace $\times (1+...) $ by $\exp (...) $ and get rid of the $z (k-1)^2$ $\endgroup$
    – Quantuple
    Jul 10, 2016 at 8:16
  • 1
    $\begingroup$ Also you should (1) calculate the time integral $\hat {V}^{(j)} $ over each path $j=1,...,M$ then (2) apply the BS formula to get $C (\hat {V}^{(j)}) $ and finally (3) average all of these to get $P_0=\frac {1}{M} \sum_{j=1}^M C (\hat {V}^{(j)})$. Istead what you did is: average all the simulated paths to compute an average path integral and apply BS. Remember that the expectation of a function is not the function of the expectation of its argument. $\endgroup$
    – Quantuple
    Jul 10, 2016 at 8:28
  • $\begingroup$ @Quantuple Thank you so much! I completely forgot about Jenson's inequality. $\endgroup$ Jul 10, 2016 at 15:35
  • $\begingroup$ @Quantuple We learnt that using the Milstein scheme produces a smaller discretisation error so I thought it'd be safer using that rather than Euler. So you propose I use the line: sigma(k) = sigma(k-1)*exp(sqrt(dt)*gammaZ(k-1) + dt*(alpha + 0.5*gammagamma)); instead? $\endgroup$ Jul 10, 2016 at 15:38
  • $\begingroup$ Your claim that Milstein scheme is better than Euler is in general true, but for the GBM case I think the difference is not significant, at least if you use Euler on log-returns. No, the equation you provide is wrong, in the last term of the exponential you should have "0.5*gamma*gamma" with a minus sign instead of a plus. Also, it's Jensen not Jenson :) $\endgroup$
    – Quantuple
    Jul 10, 2016 at 15:54

1 Answer 1


In order to compute $$ P_0 = \mathbb {E}[C (\hat{V})] $$ where $$ \hat{V} = \frac {1}{T} \int_0^T \sigma^2_s ds $$ and $$ d\sigma_t = \sigma_t (\alpha dt + \gamma dW_t) $$ using Monte Carlo, you should:

  1. Generate stochastic volatility paths over $[0,T]$ by discretising the above SDE (which here defines a GBM, not a Hull & White diffusion)
  2. Calculate the time integral $\hat {V}^{(j)} $ for each simulated path $j=1,...,M$, you can use a trapezoidal integration method for instance.
  3. Apply the BS formula to compute $C (\hat {V}^{(j)}) $
  4. Average all of the latter call prices to finally obtain $$P_0=\frac {1}{M} \sum_{j=1}^M C (\hat {V}^{(j)})$$.

Instead what you did is (see bullet points in your original question): average all simulated stochastic volatility paths to derive an "average volatility path" from which you calculate an average volatility $\hat {V} $ over $[0,T] $ and apply BS formula using the latter volatility figure.

Your approach is thus wrong. Remember that the expectation of a function of a random variable is in general not the function of the expectation of that random variable (except for linear functions obviously), see Jensen's inequality.

It is a common rookie mistake though, don't worry :)

  • $\begingroup$ is trapezoidal discretization preferable over rectangular (starting at $\sigma_0$)? slightly wondering about the implications of using $\sigma_T$ given the conditioning. $\endgroup$
    – jherek
    Jul 22, 2022 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.