function [IV, Params, RMSE, abs_residuals] = Calibration(Data, Spot, r, tau, Date, JT, underlying, expansion)
% Calibration of the Edgeworth/GramCharlier model
model = expansion;

if strcmp(underlying,'SPX')     % IC and bounds for the SPX
    if JT == 1  % Merton
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, muj, sigmaj
        init = [-0.3335, 1.0433, 0.1340, -37.3968, 0.1173, 1.0160, -0.0409, 0.0446];
        LB = [-1, 0, 0, -1000, 0, 0, -0.7, 0];
        UB = [0, 20, 0.8, 1000, 1, inf, 0.7, 0.8];
    elseif JT == 2   % Kou
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.3335, 1.0433, 0.1340, -3.3968, 0.1173, 5000, 5, 5];
        LB = [-1, 0, 0, -1000, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 1000, 1, inf, 1000, 1000];
    elseif JT == 3    % VG
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.3335, 1.0433, 0.1340, -3.3968, 0.1623, 10, 5, 5];
        LB = [-1, 0, 0, -220, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 220, 1, inf, 6000, 6000];
    else    % NIG
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.2335, 1.0433, 0.2340, -3.3968, 0.1273, 10, 5, 5];
        LB = [-1, 0, 0, -250, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 250, 1, inf, 1000, 1000];
    end
else        % IC and bounds for the DAX and Euro Stoxx 50
    if JT == 1  % Merton
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, muj, sigmaj
        init = [-0.3335, 1.0433, 0.1340, -25.3968, 0.1173, 1.0160, -0.0409, 0.0446];
        LB = [-1, 0, 0, -1000, 0, 0, -0.7, 0];
        UB = [0, 20, 0.8, 1000, 1, inf, 0.7, 0.8];
    elseif JT == 2   % Kou
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.3335, 1.0433, 0.1345, -3.3968, 0.1673, 5000, 5, 5];
        LB = [-1, 0, 0, -1000, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 1000, 1, inf, inf, inf];
    elseif JT == 3    % VG
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.3335, 1.0433, 0.1340, -37.3968, 0.2123, 10, 5, 5];
        LB = [-1, 0, 0, -450, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 450, 1, inf, inf, inf];
    else    % NIG
        % Parameteres: rho, beta, sigma, alpha, eta, lambda, psip, psim
        init = [-0.2335, 1.0433, 0.2340, -3.3968, 0.2573, 10, 5, 5];
        LB = [-1, 0, 0, -500, 0, 0, 1, 1];
        UB = [0, 20, 0.8, 500, 1, inf, inf, inf];
    end
end

options = optimoptions('lsqnonlin','FunctionTolerance',1e-10,'MaxFunctionEvaluations',5e4,'OptimalityTolerance',1e-10,'StepTolerance',1e-10,'MaxIterations',5000);
[p_exp,error_exp,residuals_exp] = lsqnonlin(@(params) objfun(Data,Spot,r,tau,params,model,JT),init,LB,UB,options);
Params.EXP = p_exp;     % Calibrated parameters
opt_number = length(Data(:,1));     % Number of options
RMSE.EXP = 100*sqrt(error_exp/opt_number);      % Root mean squared error

% Edgeworth/GramCharlier IV
IV.EXP =  objfun([Data(:,1:end-1), zeros(opt_number,1)],Spot,r,tau,p_exp,model,JT); 

% Calibration of the BSM model
model = 'BSM';

if JT == 1      % Merton
    % Parameteres: sigma, lambda, muj, sigmaj which are in positions 3,6,7,8
    init = [0, 0, 0.2340, 0, 0, 1.0160 -0.0409 0.0446];
    LB = [0, 0, 0, 0, 0, 0, -0.5, 0];
    UB = [0, 0, 0.8, 0, 0, inf, 0.5, 0.8];
else    % Tempered stable jumps
    % Parameteres: sigma, lambda, psip, psim which are in positions 3,6,7,8
    init = [0, 0, 0.2340, 0, 0, 1.0160, 20, 20];
    LB = [0, 0, 0, 0, 0, 0, 1, 1];
    UB = [0, 0, 0.8, 0, 0, inf, inf, inf];
end

[p_bsm, error_bsm, residuals_bsm] = lsqnonlin(@(params) objfun(Data,Spot,r,tau,params,model,JT),init,LB,UB,options);
Params.BSM = p_bsm([3,6:8]);    % Calibrated parameters
RMSE.BSM = 100*sqrt(error_bsm/opt_number);      % Root mean squared error

% Absolute residuals of the calibrations
abs_residuals = 100*[abs(residuals_bsm), abs(residuals_exp)];

% BSM IV
IV.BSM =  objfun([Data(:,1:end-1), zeros(opt_number,1)],Spot,r,tau,p_bsm,model,JT); 

% Plot of Market, Edgeworth/GramCharlier and BSM IV
figure()
plot(Data(:,4),Data(:,5),'.','Markersize', 22)
hold on
grid on
plot(Data(:,4),IV.EXP,'LineWidth',1.2)
plot(Data(:,4),IV.BSM,'--','LineWidth',1.2)
title(strcat('Option pricing fit on'," ", string(Date)))
xlabel('Moneyness')
ylabel('Implied Volatility')
legend('Market',expansion,'BSM')
set(findall(gca,'Type','Line'),'Linewidth',2)
set(gca,'Fontsize',16)

end