%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Driver for parameter and implied volatilty fitting under Merton jump-diffusion
%%% model
%%%
%%% Copyright (C) 2008  Vladimir Surkov
%%%
%%%
%%% This program is free software: you can redistribute it and/or modify
%%% it under the terms of the GNU General Public License as published by
%%% the Free Software Foundation, either version 3 of the License, or
%%% (at your option) any later version.
%%%
%%% This program is distributed in the hope that it will be useful,
%%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%% GNU General Public License for more details.
%%% 
%%% You should have received a copy of the GNU General Public License
%%% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [sigma, lambda, mu, nu] = Driver_ParameterFitting_MertonJumpDiffusion()

load('SPY_Option_Prices');
K = SPY_Option_Prices(:,1)';
C = SPY_Option_Prices(:,2)';
P = SPY_Option_Prices(:,3)';

S = 117.5;
T = 1/12;
r = 0.02;
q = 0.0221;


f = @(x)(PricingError(S, K, T, x(1), x(2), x(3), x(4), r, q, C, P));
x_l = [0.1, 0.1, -0.5, 0.01];
x_0 = [0.2, 1.0, -0.1, 0.1];
x_u = [0.5, 5, 0.5, 0.3];

x = lsqnonlin(f, x_0, x_l, x_u);
sigma = x(1); lambda = x(2); mu = x(3); nu = x(4);

[values_call, values_put] = CarrMadanFFT_European_MertonJumpDiffusion(S, K, T, sigma, lambda, mu, nu, r, q);
[sigma_call_model, sigma_put_model] = ImpliedVolatility(S, K, T, r, q, values_call, values_put);
[sigma_call_market, sigma_put_market] = ImpliedVolatility(S, K, T, r, q, C, P);

subplot(2,1,1);
plot(K, sigma_call_model, 'r', K, sigma_call_market, 'b');
subplot(2,1,2);
plot(K, sigma_put_model, 'r', K, sigma_put_market, 'b');

end

function error = PricingError(S, K, T, sigma, lambda, mu, nu, r, q, C, P)
    w = blsgamma(S, K, r, T, sigma, q);
    [values_call, values_put] = CarrMadanFFT_European_MertonJumpDiffusion(S, K, T, sigma, lambda, mu, nu, r, q);
    error = sqrt([(values_call - C)./w (values_put - P)./w]);
end