%---1) MONTECARLO: SIMULIAMO nTrials TRAIETTORIE DI PREZZO in HESTON e in
% BLACK & SCHOLES ---%%% Prezzo, volatilità, curtosi

clc;
clear all;
warning off;
close all;

nTrials = 200;
strike = 15;
inizioPrezzi = 5;
numPrezzi = 1;

speed = 3;
nPeriods = 500;
scadenza = 50;
mu = 0.02;
deltaTime=1/nPeriods;
volavola = 3;
level = 0.9;
startstatevola = 1.3;
call = 1;
r = 0.01;
correlation = 0.5;

startState = [15; startstatevola];
% B&S: usa vola di lp di Heston (con X)
[PrezziBeS, VolaBeS, TempiBeS] = generaPrezziBeS(mu, level, startState, nPeriods, nTrials, deltaTime);
% Heston
[PrezziHeston, VolaHeston, tempiHeston] = generaPrezziHeston(mu, speed, level, volavola, correlation, startState, nPeriods, nTrials, deltaTime);
volaHeston = mean(VolaHeston, 1); % vettore con le volatilità medie nelle diverse scadenze

% Plotta prezzi e volatilità dei due
figure()
subplot(2,2,1);
hold on
for i = 1:nTrials
plot(TempiBeS, PrezziBeS(i,:));%prezzo
title('B&S: PREZZI');
end
subplot(2,2,3);
hold on;
for i = 2:nTrials
plot(TempiBeS, VolaBeS(i,:)./PrezziBeS(i,:));%prezzo
title('B&S: sigma');
end
subplot(2,2,2);
hold on
for i = 1:nTrials
plot(tempiHeston, PrezziHeston(i,:));%prezzo
title('HESTON: PREZZI');
end
subplot(2,2,4);
hold on;
for i = 2:nTrials
plot(tempiHeston, VolaHeston(i,:));
title('Heston: sigma');
end
% MOSTRA LA CURTOSI
figure();
lineWidth = 2;
subplot(1,2,1);
hold on;
returns = tick2ret(PrezziBeS'); % Calcola rendimenti delle traiettorie generate con Heston
muHist = mean(mean(returns));
sigmaHist = sqrt(var(returns(:)));
curtosi = kurtosis(returns(:));
histogram(real(returns),'Normalization','pdf');
y = -0.15:0.001:0.15;
f = exp(-(y-muHist).^2./(2*sigmaHist^2))./(sigmaHist*sqrt(2*pi));
plot(y,f,'LineWidth',lineWidth) % Disegna la normale
title(strcat('Black & Scholes, curtosi: ', num2str(curtosi)));

subplot(1,2,2);
hold on;
returns = tick2ret(real(PrezziHeston'));
muHist = mean(mean(returns));
sigmaHist = sqrt(var(returns(:)));
curtosi = kurtosis(returns(:));
histogram(real(returns),'Normalization','pdf');
y = -0.15:0.001:0.15;
f = exp(-(y-muHist).^2./(2*sigmaHist^2))./(sigmaHist*sqrt(2*pi));
plot(y,f,'LineWidth',lineWidth)
title(strcat('Heston, curtosi: ', num2str(curtosi)));

%% 2) MONTECARLO: SMILE (prezzo fisso-strike varia) HESTON

clc;
clear all;
% close all;
rng default;

%parametri da modificare
nTrials = 400
nSteps = 100;
prezzo = 15;
inizioStrike = 10;
numStrike = 10;

%parametri da lasciare
speed = 0.2;
nPeriods = 5;
deltaTime=1/nPeriods;

ignoraVola =true;

mu = 0.02;
volavola = 1.83;
levola = 0.02; %volatilità di partenza e di lp
level = levola;
startstatevola = levola;

r = 0.01; % tasso istantaneo di rendimento risk free
correlation = -0.5;

Vola3Scadenze = zeros(numStrike,3);
scadenza1 = 2;
scadenza2 = 3;
scadenza3 = 4;
prezziScad = zeros(numStrike,1);
rng default;
startState = [prezzo; startstatevola];

for u = 1:numStrike
[PrezziHeston, tempiHeston] = generaPrezziHeston2(mu, speed, level, volavola, correlation, startState, nPeriods, nTrials, deltaTime, ignoraVola, nSteps);

if (inizioStrike+u)>prezzo
call = 1;
callbls = true;
else
call = 0;
callbls = false;
end

% interrompi(u)
payoffsHeston = optionPrices2(call, inizioStrike+u, r, 1/nPeriods, tempiHeston, PrezziHeston);
mediePayoffs = payoffsHeston;

volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza1/nPeriods, mediePayoffs(scadenza1),200, 0,[],callbls);
vola3Scadenze(u,1) = volatility;
volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza2/nPeriods, mediePayoffs(scadenza2),200, 0,[],callbls);
vola3Scadenze(u,2) = volatility;
prezziScad(u) = mediePayoffs(scadenza2);
volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza3/nPeriods, mediePayoffs(scadenza3),200, 0,[],callbls);
vola3Scadenze(u,3) = volatility;
end
% beep on;
% beep;
figure();
scatter(inizioStrike+1:numStrike+inizioStrike, prezziScad);
title(strcat('Prezzi opzioni put - call, HESTON, So= ', num2str(prezzo), 'maturity = ', num2str(scadenza2/nPeriods*250), ' giorni'));
xlabel('Strike');
ylabel('Prezzo');

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,1),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,1),'b','x');
a = vola3Scadenze(1:numStrike,1)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),2);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, HESTON, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza1/nPeriods*250), ' giorni'));
ylabel('Implied volatility');
xlabel('Strike');
legend(strcat('Y=(', num2str(p(1)),')X^2 +(', num2str(p(2)),')X+(', num2str(p(3)),')'));

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,2),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,2),'b','x');
a = vola3Scadenze(1:numStrike,2)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),2);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, HESTON, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza2/nPeriods*250), ' giorni'));
ylabel('Implied volatility');
xlabel('Strike');
legend(strcat('Y=(', num2str(p(1)),')X^2 +(', num2str(p(2)),')X+(', num2str(p(3)),')'));

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,3),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,3),'b','x');
a = vola3Scadenze(1:numStrike,3)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),2);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, HESTON, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza3/nPeriods*250), ' giorni'));
ylabel('Implied volatility');
xlabel('Strike');
legend(strcat('Y=(', num2str(p(1)),')X^2 +(', num2str(p(2)),')X+(', num2str(p(3)),')'));

%% MONTECARLO: SMILE (B&S)

clc;
clear all;
% close all;
rng default;

%parametri da modificare
nTrials = 1000;

prezzo = 15;
inizioStrike = 10;
numStrike = 10;

%parametri da lasciare
speed = 0.2;
nPeriods = 5;
deltaTime=1/nPeriods;
nSteps = 100;
ignoraVola =true;

mu = 0.02;
volavola = 1.83;
levola = 0.5;
level = levola;
startstatevola = levola;
r = 0.01;
correlation = -0.5;

Vola3Scadenze = zeros(numStrike,3);
scadenza1 = 2;
scadenza2 = 3;
scadenza3 = 4;
prezziScad = zeros(numStrike,1);
rng default;
startState = [prezzo; startstatevola];

for u = 1:numStrike
[PrezziBeS, VolaBeS, tempiBeS] = generaPrezziBeS2(mu, level, startState, nPeriods, nTrials, deltaTime, ignoraVola, nSteps);
clear VolaBeS;
if (inizioStrike+u)>prezzo
call = 1;
callbls = true;
else
call = 0;
callbls = false;
end
u
payoffsBeS = optionPrices(call, inizioStrike+u, r, 1/nPeriods, tempiBeS, PrezziBeS);
mediePayoffs = payoffsBeS;

volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza1/nPeriods, mediePayoffs(scadenza1),200, 0,[],callbls);
vola3Scadenze(u,1) = volatility;

volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza2/nPeriods, mediePayoffs(scadenza2),200, 0,[],callbls);
vola3Scadenze(u,2) = volatility;
prezziScad(u) = mediePayoffs(scadenza2);

volatility = blsimpv(prezzo, inizioStrike+u, r, scadenza3/nPeriods, mediePayoffs(scadenza3),200, 0,[],callbls);
vola3Scadenze(u,3) = volatility;
end
% beep on;
% beep;
figure();
scatter(inizioStrike+1:numStrike+inizioStrike, prezziScad);
title(strcat('Prezzi opzioni put - call, B&S, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza2/nPeriods*250), ' giorni'));
xlabel('Strike');
ylabel('Prezzo');

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,1),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,1),'b','x');
a = vola3Scadenze(1:numStrike,1)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),1);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, B&S, So= ', num2str(prezzo), 'maturity= ', num2str(scadenza1/nPeriods*250), ' giorni'));
xlabel('Implied volatility');
ylabel('Strike');

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,2),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,2),'b','x');
a = vola3Scadenze(1:numStrike,2)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),1);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, B&S, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza2/nPeriods*250), ' giorni'));
xlabel('Implied volatility');
ylabel('Strike');

figure();
hold on;
scatter((inizioStrike+1:1:prezzo), vola3Scadenze(1:prezzo-inizioStrike,3),'r','x');
scatter((prezzo+1:1:numStrike+inizioStrike), vola3Scadenze(prezzo-inizioStrike+1:numStrike,3),'b','x');
a = vola3Scadenze(1:numStrike,3)';
ind = inizioStrike+1:numStrike+inizioStrike;
k = ~isnan(a);
p = polyfit(ind(k),a(k),1);
refcurve(p);
axis([10 20 0 0.8]);
title(strcat('Moneyness - implied volatility, B&S, So= ', num2str(prezzo), ' maturity= ', num2str(scadenza3/nPeriods*250), ' giorni'));
xlabel('Implied volatility');
ylabel('Strike');
% beep;

%% 3) - FORMULA ANALITICA - volatility surface
clear all
% close all;
s0=1500;
vola0 = 1;
level = 0.1;
volavola = 5;
mu = 0.05;
speed = 2;
correlation = -0.3;
r = 0.03;

inizioK = 500;
fineK = 2500;
nintervalliK = 40;
intervalliK = ceil(fineK-inizioK)/nintervalliK;

k = inizioK:intervalliK:fineK;
k = k/s0;

ttm = 2:-0.04:0.1;

[X, Y] = meshgrid(ttm, k);
Z = zeros(length(k),length(ttm));
for i=1:length(k)
for j=1:length(ttm)
OptionPrice = hestonAnalitica(s0/s0, vola0, speed, level, volavola,...
correlation, r, ttm(j), k(i));
Z(i, j) = blsimpv(s0/s0, k(i), r, ttm(j), max(OptionPrice,0),100);
end
i
end

Y=-Y;
X=X*250;
figure();
mesh(X, Y, Z);
xlabel('Time to Maturity: giorni');
ylabel('Moneyness: Strike/prezzo');
zlabel('Implied volatility');
title('VOLATILITY SURFACE');

%% 4) CALIBRAZIONE
%shuffledArray = orderedArray(randperm(size(orderedArray,1)),:);
clear all; close all;
%variabili da prendere nella funzione di costo
global strikes; global maturity; global price; global r;

data = load('c_data2.csv');
% data = data(randperm(size(data,1)),:);
% for u = 5:-1:0
u = 0;
numData1 = u*380+1;
numData2 = numData1+30;

% numData1 = 1;
% numData2 = 30;
% numData = length(data(:,1));
strikes = data(numData1:numData2,1);
maturity = data(numData1:numData2,2); %frazioni di anno fra 3/9/2017 e scadenza
bid = data(numData1:numData2,3);
ask = data(numData1:numData2,4);
price = data(numData1:numData2,5);
% x0:
% s0, x1:vola0, x2:speed, x3:level, x4:volavola,x5:correlation, x6:r, ttm, k
% x0 = [1;0.5;0.5;0.5;-0.5; 0.01]
x0 = [0.5;0.5;0.5;0.5;-0.5]
r = 0.01;
startprice = 2476.55;

%price = hestonAnalitica(startprice, x0(1),x0(2),x0(3),x0(4),x0(5),x0(6), ttm, k)

%c = costFunction(x0)

lb = [0;0;0;0;-1]
ub = [1;1;10;10;1]
x = lsqnonlin(@costFunction,x0,lb,ub)

n = length(price)
predicted_price = zeros(n,1);

for i=1:n
predicted_price(i) = hestonAnalitica(startprice, x(1),x(2),x(3),x(4),x(5),r, maturity(i), strikes(i));
end
%[bid,price, predicted_price,ask]
subplot(3,2,u+1);
scatter(price,predicted_price,'filled');
hold on;
% confronto predetti/bisettrice
plot([min(price),max(price)],[min(price),max(price)]);
title(strcat(num2str(min(strikes)),' - ',num2str(max(strikes))))
xlabel('predetti')
ylabel('effettivi')
%%per far respirare il pc
% if mod(u,2)==1
% beep();
% disp('premere un tasto');
% pause();
% end;
% end

%%

inizioK = 1500;
fineK = 3000;
nintervalliK = 40;
intervalliK = ceil(fineK-inizioK)/nintervalliK;
s0 = startprice;

k = inizioK:intervalliK:fineK;
k = k/s0;

ttm = 1:-0.02:0.05;

[X, Y] = meshgrid(ttm, k);
Z = zeros(length(k),length(ttm));
for i=1:length(k)
for j=1:length(ttm)
OptionPrice = hestonAnalitica(s0/s0, x(1), x(2), x(3), x(4), x(5), r, ttm(j), k(i));
Z(i, j) = blsimpv(s0/s0, k(i), r, ttm(j), max(OptionPrice,0),100);
end
i
end

Y=-Y;
X=X*250;
figure();
mesh(X, Y, Z);
xlabel('Time to Maturity: giorni');
ylabel('Moneyness: Strike/prezzo');
zlabel('Implied volatility');
title('VOLATILITY SURFACE');

function [pay] = optionPrices(call, strike, r, deltaTime, tempi, prezziGenerati)
%HESTONOPTIONPRICE Summary of this function goes here

attualizzazione = exp(-tempi*r);
payoffs = prezziGenerati-strike;
if call==1
payoffs = max(payoffs, 0);
else
payoffs = max (-payoffs, 0);
end
% Attualizza i payoff
pay = mean(payoffs);
pay = pay*diag(attualizzazione);

end

function [ prezzo ] = hestonAnalitica(s0, vola0, speed, level, volavola,...
correlation, r, ttm, k )
%HESTONANALITICA Summary of this function goes here
% Detailed explanation goes here

integrando1 = @(s,s0, vola0, level, speed, volavola, r, correlation, ttm, k) real(exp(-i.*s*log(k)).*...
funzCarHeston(s0, vola0, level, speed, volavola, r, correlation, ttm, s-i)./...
(i*s.*funzCarHeston(s0, vola0, level, speed, volavola, r, correlation, ttm, -i)));
integrando2 = @(s,s0, vola0, level, speed, volavola, r, correlation, ttm, k) real(exp(-i.*s*log(k)).*...
funzCarHeston(s0, vola0, level, speed, volavola, r, correlation, ttm, s)./...
(i*s));

integrale1 = integral(@(s)integrando1(s,s0, vola0, level, speed, volavola, r, correlation, ttm, k),0,100);
integrale2 = real(integral(@(s)integrando2(s,s0, vola0, level, speed, volavola, r, correlation, ttm, k),0,100));

p1 = 0.5+integrale1/pi;
p2 = 0.5+integrale2/pi;

prezzo = s0*p1-k*p2*exp(-r*ttm);
end

function [prezziGenerati, volaGenerati, tempi] = generaPrezziBeS(mu, level, startState, nPeriods, nTrials, deltaTime)

% GENERAPREZZIB&S Summary of this function goes here
% Genera nTrials traiettorie secondo B&S, su nPeriods scadenze.

F = @(t, X) mu*X;
G = @(t, X) X*level;
besSde = sde(F, G, 'StartState', startState(1));
[traiettorie, tempi, z] = simByEuler(besSde, nPeriods, 'nTrials', nTrials, 'DeltaTime', deltaTime);
traiettorie = permute(traiettorie, [3, 1, 2]);
prezziGenerati = traiettorie(:, :, 1);
volaGenerati = G(tempi, prezziGenerati);

end

function costo = costFunction(x)

global strikes; global maturity; global price; global r;
%COSTFUNCTION Summary of this function goes here
% Detailed explanation goes here

startprice = 2476.55;
% k = 1925;
% ttm = 0.0383561643835616;
% actualprice = 551.15;

n = length(strikes);
cost = zeros(n,1);
for i=1:n
prezzoHest = hestonAnalitica(startprice, x(1),x(2),x(3),x(4),x(5),r, maturity(i), strikes(i));
%prezzoHest = hestonAnalitica(startprice, x(1),x(2),(x(4)+x(3)^2)/(2*x(2)),x(4),x(5),x(6), maturity(i), strikes(i));
cost(i) = (prezzoHest-price(i)).^2;
end
% costo = sum(cost)/n;
costo = cost;
end

function [ funzCarHeston ] = funzCarHeston(s0, vola0, level, speed, volavola, r, correlation, ttm, s)
%FUNZCARHESTON Summary of this function goes here
% Detailed explanation goes here

alpha = -s.*s/2 - i*s/2;
beta = alpha - correlation*volavola*i*s;
gamma = volavola*volavola/2;
h = sqrt(beta.*beta - 4*alpha*gamma);
rplus = (beta + h)/volavola/volavola;
rminus = (beta - h)/volavola/volavola;
g=rminus./rplus;

C = spee

function [prezziGenerati, tempi] = generaPrezziHeston(mu, speed, level, volavola, correlation, startState, nPeriods, nTrials, deltaTime, ignoraVola, nSteps)
% GENERAPREZZIHESTON Summary of this function goes here
% Genera nTrials traiettorie secondo il processo di Heston, su nPeriods scadenze.

hestonSde = heston(mu, speed, level, volavola,...
'StartState', startState, 'Correlation', correlation);
[traiettorie, tempi] = simByEuler(hestonSde, nPeriods,...
'nTrials', nTrials, 'DeltaTime', deltaTime, 'nSteps', nSteps);
prezziGenerati = squeeze(traiettorie(:,1,:))';
end

d*(rminus*ttm-(2/volavola^2).*log((1-g.*exp(-h*ttm))./(1-g)));
D = rminus.*(1-exp(-h*ttm))./(1-g.*exp(-h*ttm));

funzCarHeston = exp(C*level+D*vola0+i*s*log(s0*exp(r*ttm)));

end