# Hull White Stochastic Volatility Model in Matlab

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

• 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$ – Quantuple Jul 10 '16 at 8:16
• 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. – Quantuple Jul 10 '16 at 8:28
• @Quantuple Thank you so much! I completely forgot about Jenson's inequality. – Dabshffabjvs Jul 10 '16 at 15:35
• @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? – Dabshffabjvs Jul 10 '16 at 15:38
• 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 :) – Quantuple Jul 10 '16 at 15:54

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.