Skip to content

Commit

Permalink
addes files, first commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@skiamu
skiamu committed Aug 20, 2017
0 parents commit b7baa9b
Show file tree
Hide file tree
Showing 37 changed files with 2,974 additions and 0 deletions.
114 changes: 114 additions & 0 deletions Call_Operator_Splitting_Levy_VG.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
function price=Call_Operator_Splitting_Levy_VG
%% Call_Operator_Splitting_Levy
% European Call, VG process

% Model parameters
T = 1;S0 = 30;r = 0.1;K = 30;
nu=0.16; etaN=0.0694; etaP=0.0166; %Levy density parameters
% param = [6.25;60.2;14.4];
% Discretization parameters
N = 200; M = 100;
% Grids
dt=T/M;
Smin=10; Smax=300;
xmin=log(Smin/K); xmax=log(Smax/K);
x=linspace(xmin,xmax,N+1); dx=x(2)-x(1);
nodes=x(2:end-1)'; % interior nodes

% Compute sigma_Eps, alpha, lambda
epsilon=0.01;
k=@(y) 1./(nu*abs(y)).*exp(-abs(y)/etaN).*(y<0)+...
1./(nu*abs(y)).*exp(-abs(y)/etaP).*(y>0);
% model = 'VG2';
% k = LevyDensity(param, model);
[alpha,lambda,sigma,Bl,Br]=Levy_Integrals(k,N,epsilon)

% Matrix
Dd=-(sigma^2/(2*(dx^2))+(sigma^2/2+alpha)/(2*dx));
Dm=1/dt-(-sigma^2/(dx^2)-lambda);
Du=-(sigma^2/(2*(dx^2))-(sigma^2/2+alpha)/(2*dx));
e=ones(N-1,1);
M1=spdiags([Dd*e Dm*e Du*e],-1:1,N-1,N-1);

% Forward-in-time loop
u=max( exp(nodes)-1, 0); % initial solution
BC=zeros(N-1,1);
for j=1:M
% Boundary condition
BC(end)=-Du*(exp(xmax)-1);
% Compute I
I=Levy_Int(k,nodes,u,Bl,Br,N,epsilon);
%Solve the linear system
u=M1\(u/dt+BC+I);
end
% Plot and interpolation
S=K*exp(nodes-r*T);
C=K*u*exp(-r*T);
figure
plot(S,C);
price=interp1(S,C,S0,'spline')
if lambda==0
price_ex=blsprice(S0,K,r,T,sigma)
end
end

function [alpha,lambda,sigma,Bl,Br]=Levy_Integrals(k,N,epsilon)
Tol=10^-12; Bl=-0.01; Br=0.01;
% while (k(Bl)>Tol)
% Bl=Bl-0.002;
% end
% while (k(Br)>Tol)
% Br=Br+0.002;
% end
Bl = fsolve(@(y) k(y) - Tol, Bl);
Br = fsolve(@(y) k(y) - Tol, Br);
qnodes=linspace(-epsilon,epsilon,2*N);
sigma=trapz(qnodes,(qnodes.^2).*k(qnodes));
qnodes1=linspace(Bl,-epsilon,N);
qnodes2=linspace(epsilon,Br,N);
lambda=trapz(qnodes1,k(qnodes1))+trapz(qnodes2,k(qnodes2));
alpha=trapz(qnodes1,(exp(qnodes1)-1).*k(qnodes1))+...
trapz(qnodes2,(exp(qnodes2)-1).*k(qnodes2));
figure; hold on
plot(qnodes1,k(qnodes1),'b');
plot(qnodes2,k(qnodes2),'b');
plot(qnodes,k(qnodes),'r');


end

function I=Levy_Int(k,nodes,u,Bl,Br,N,epsilon)
qnodes1=linspace(Bl,-epsilon,round(N/2));
dq1=qnodes1(2)-qnodes1(1);
qnodes2=linspace(epsilon,Br,round(N/2));
dq2=qnodes2(2)-qnodes2(1);
I=zeros(N-1,1);
for i=1:N-1
I(i)=trapz(f_u(nodes(i)+qnodes1,nodes,u).*k(qnodes1))*dq1+...
trapz(f_u(nodes(i)+qnodes2,nodes,u).*k(qnodes2))*dq2;
end
end

function f=f_u(z,nodes,u)
f=zeros(size(z));
index=find(z>nodes(end));
f(index)=exp(z(index))-1;
index=find( (z>=nodes(1)).*(z<=nodes(end)) );
f(index)=interp1( nodes,u,z(index));












end




112 changes: 112 additions & 0 deletions Call_Operator_Splitting_Levy_VG.m~
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
function price=Call_Operator_Splitting_Levy_VG
%% Call_Operator_Splitting_Levy
% European Call, VG process

% Model parameters
K=100; S0=100; r=0.01; T=1;
% nu=0.16; etaN=0.0694; etaP=0.0166; %Levy density parameters
param = [6.25;60.2;14.4];
% Discretization parameters
N=2000; M=100;
% Grids
dt=T/M;
Smin=10; Smax=300;
xmin=log(Smin/K); xmax=log(Smax/K);
x=linspace(xmin,xmax,N+1); dx=x(2)-x(1);
nodes=x(2:end-1)'; % interior nodes

% Compute sigma_Eps, alpha, lambda
epsilon=0.01;
% k=@(y) 1./(nu*abs(y)).*exp(-abs(y)/etaN).*(y<0)+...
% 1./(nu*abs(y)).*exp(-abs(y)/etaP).*(y>0);
model =
k = LevyDensity(param, model);
[alpha,lambda,sigma,Bl,Br]=Levy_Integrals(k,N,epsilon);

% Matrix
Dd=-(sigma^2/(2*(dx^2))+(sigma^2/2+alpha)/(2*dx));
Dm=1/dt-(-sigma^2/(dx^2)-lambda);
Du=-(sigma^2/(2*(dx^2))-(sigma^2/2+alpha)/(2*dx));
e=ones(N-1,1);
M1=spdiags([Dd*e Dm*e Du*e],-1:1,N-1,N-1);

% Forward-in-time loop
u=max( exp(nodes)-1, 0); % initial solution
BC=zeros(N-1,1);
for j=1:M
% Boundary condition
BC(end)=-Du*(exp(xmax)-1);
% Compute I
I=Levy_Int(k,nodes,u,Bl,Br,N,epsilon);
%Solve the linear system
u=M1\(u/dt+BC+I);
end
% Plot and interpolation
S=K*exp(nodes-r*T);
C=K*u*exp(-r*T);
figure
plot(S,C);
price=interp1(S,C,S0,'spline')
if lambda==0
price_ex=blsprice(S0,K,r,T,sigma)
end
end

function [alpha,lambda,sigma,Bl,Br]=Levy_Integrals(k,N,epsilon)
Tol=10^-10; Bl=-1; Br=1;
while (k(Bl)>Tol)
Bl=Bl-2;
end
while (k(Br)>Tol)
Br=Br+2;
end
qnodes=linspace(-epsilon,epsilon,2*N);
sigma=trapz(qnodes,(qnodes.^2).*k(qnodes));
qnodes1=linspace(Bl,-epsilon,N);
qnodes2=linspace(epsilon,Br,N);
lambda=trapz(qnodes1,k(qnodes1))+trapz(qnodes2,k(qnodes2));
alpha=trapz(qnodes1,(exp(qnodes1)-1).*k(qnodes1))+...
trapz(qnodes2,(exp(qnodes2)-1).*k(qnodes2));
figure; hold on
plot(qnodes1,k(qnodes1),'b');
plot(qnodes2,k(qnodes2),'b');
plot(qnodes,k(qnodes),'r');


end

function I=Levy_Int(k,nodes,u,Bl,Br,N,epsilon)
qnodes1=linspace(Bl,-epsilon,round(N/2));
dq1=qnodes1(2)-qnodes1(1);
qnodes2=linspace(epsilon,Br,round(N/2));
dq2=qnodes2(2)-qnodes2(1);
I=zeros(N-1,1);
for i=1:N-1
I(i)=trapz(f_u(nodes(i)+qnodes1,nodes,u).*k(qnodes1))*dq1+...
trapz(f_u(nodes(i)+qnodes2,nodes,u).*k(qnodes2))*dq2;
end
end

function f=f_u(z,nodes,u)
f=zeros(size(z));
index=find(z>nodes(end));
f(index)=exp(z(index))-1;
index=find( (z>=nodes(1)).*(z<=nodes(end)) );
f(index)=interp1( nodes,u,z(index));












end




131 changes: 131 additions & 0 deletions CharacteristicExp.m
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
function [ Psi ] = CharacteristicExp( model, param )
%Characteristic_Exp restituisce l'esponente caratteristico del modello
%indicato in input
%
% INPUT:
% model = striga che indica il modello di cui si vuole l'espoinente
% caratteristico
%
% param = vettore parametri del modello
% --> 'Heston'
% param.k;
% param.theta;
% param.epsilon;
% param.rho;
% param.V0;
%
% --> 'BS'
% param(1) = sigma
% --> 'Kou'
% param.sigma;
% param.lambda;
% param.lambda_p;
% param.lambda_m;
% param.p;
%
% --> 'Merton'
% param(1) = lambda;
% param(2) = sigma;
% param(3) = mu;
% param(4) = delta;
%
% --> 'VG'
% param(1) = sigma
% param(2) = kappa
% OUTPUT:
% Psi = esponente caratteristico (function handle)


switch model

case 'Heston'

% parametri
k = param.k;
theta = param.theta;
rho = param.rho;
epsilon = param.epsilon;
V0 = param.V0;

gamma = @(u) k - rho * epsilon .* u * 1i;

alpha = @(u) - 0.5 .* (u.^2 + 1i .* u);

d = @(u) sqrt( gamma(u).^2 - 2 * (epsilon^2) * alpha(u));

C = @(t,u) (-2 * k * theta / (epsilon^2)) .* ...
(2 * log( (2 * d(u) - (d(u) - gamma(u)) .* ...
(1 - exp(-d(u) * t)))./ (2*d(u)) ) + (d(u) - gamma(u)) * t );

B = @(t,u) ((2 * alpha(u) .* (1 - exp(- d(u) * t))) * V0) ./ ...
((2 * d(u) - (d(u) - gamma(u)) .* (1 - exp(-d(u) * t))));

Psi = @(t,u) B(t,u) + C(t,u);

case 'Merton'
% estrazione parametri
% parametri Merton: 1) sigma = diffusion volatility
% 2) lambda = jump intensity
% 3) mu = mean jump
% 4) delta = jump std
sigma = param(1);
lambda = param(2);
mu = param(3);
delta = param(4);

% è l'esponente caratteristico con drift nullo
Psi = @(u) - 0.5 * u.^2 * sigma^2 + lambda * (exp(- 0.5 * delta^2 *...
u.^2 + 1i * mu * u) - 1);

case 'Kou'
% estrazione parametri
% parametri Kou: 1) sigma = diffusion volatility;
% lambda = jump intensity
% lambda_p = positive jump parameter
% param_m = negative jump parameter
% p = probability positive jump
sigma = param(1);
lambda = param(2);
lambda_p = param(3);
lambda_m = param(4);
p = param(5);

% esponente caratteristico con drift nullo
Psi = @(u) (-sigma^2 .* u.^2) / 2 + 1i .* u .* lambda .* (p ./ (lambda_p - ...
1i .* u) - (1 - p) ./ (lambda_m + 1i .* u));

case 'NIG'
sigma = param(1);
kappa = param(2);

% esponente caratteristico con drift nullo
Psi = @(u) 1 / kappa - (1 / kappa) * sqrt( 1 + u.^2 * sigma^2 * kappa);

case 'VG'
% parametri VG: 1) sigma = diffusion volatility
% 2) kappa = subordinator volatility
sigma = param(1);
kappa = param(2);

% esponente caratteristico con drift nullo
Psi = @(u) - (1 / kappa) * log(1 + 0.5 * u.^2 * sigma^2 * kappa);

case 'BS'
sigma = param(1);

% esponente caratteristico con drift nullo
Psi = @(v) - (sigma^2) / 2 .* v.^2;
case 'alpha'
alpha = param(1);
sigma = param(2);

% esponente caratteristico con drift nullo
Psi = @(u) -sigma^alpha * abs(u).^alpha;
otherwise
error('model invalid : %s', model);

end


end

Loading

0 comments on commit b7baa9b

Please sign in to comment.