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:
- Simulate Paths of sigma using the Milstein scheme for discretising an SDE
- Average points at each time incriment for all the simulations and compute \bar{V}.
- Find the Black Scholes Price
- 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