Hull White Stochastic Volatility Model in Matlab

Question

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).

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

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

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)));
     end
     sigma_path(:,j)=sigma;
end
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;
end

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 :)