|
|
 |
|
|
Financial product design and pricing |
|
|
經歷
興趣
Top |
|
Top |
|
|
Top |
|
Top |
|
HW1
function X=BSCall(S,K,r,sigma,T)
Q=(log(S/K)+(r+sigma^2/2)*T)/(sigma*sqrt(T));
W=Q-(sigma*sqrt(T));
X=S*normcdf(Q)-K*(exp(-r*T)*normcdf(W));
function X=BSPut(S,K,r,sigma,T)
Q=(log(S/K)+(r+sigma^2/2)*T)/(sigma*sqrt(T));
W=Q-(sigma*sqrt(T));
X=-S*normcdf(-Q)+K*(exp(-r*T)*normcdf(-W));
function B = PV(y,n,C)
for i=1:n
D(i)= 1/((1+y)^(i-1));
end
D*C
HW2 Deltax=50; r=0.05; t=0.5:0.5:10; sig=0.2; div=0; s=60; call_delta1= blsdelta(s,x,r,t,sig,div); call_delta2= blsdelta(s-10,x,r,t,sig,div); call_delta3= blsdelta(s-20,x,r,t,sig,div); plot(t,call_delta1,'-') hold on; plot(t,call_delta2,'-r') plot(t,call_delta3,'-g') title('Delta vs. T') xlabel('T') ylabel('Delta') legend('in the money','at the money','out of the money') axis([0 10 0 1.2]); Gamma
x=50; r=0.05; t=0.5; sig=0.2; div=0; s=30:120; call_gamma = blsgamma(s,x,r,t,sig,div); plot(s,call_gamma,'-') title('Gamma vs. s') xlabel('s') ylabel('Gamma') legend('Gamma') Theta
x=50; r=0.05; t=0.5; sig=0.2; div=0; s=30:120; [call_theta,put_theta]=blstheta(s,x,r,t,sig,div); plot(s,call_theta,'-') title('Theta vs. s') xlabel('s') ylabel('Theta') legend('Theta') Vega
x=50; r=0.05; t=0.5; sig=0.2; div=0; s=30:120; call_vega = blsvega(s,x,r,t,sig,div); plot(s,call_vega,'-') title('Vega vs. s') xlabel('s') ylabel('Vega') legend('Vega')
HW4
%Callprice.m
function [C,P] = blsprice(S0,X,r,T,sigma)
X=20:1:80;
S0=50;
r=0.08;
sigma=0.4;
T=1;
[C,P] = blsprice(S0,X,r,T,sigma);
plot(X,C,'-r');
hold on;
T=0;
[C,P] = blsprice(S0,X,r,T,sigma);
plot(X,C,'-blue');
hold on;
axis([20 80 0 25]);
title('Call price v.s. Exercise Price');
xlabel('Exercise price');
ylabel('Call price');
legend('EurCall value','內含價值(T=0)');
function [C,P] = blsprice(S0,X,r,T,sigma)
X=20:1:80;
S0=50;
r=0.08;
sigma=0.4;
T=1;
[C,P] = blsprice(S0,X,r,T,sigma);
plot(X,P,'-r');
hold on;
T=0;
[C,P] = blsprice(S0,X,r,T,sigma);
plot(X,P,'-blue');
hold on;
axis([20 80 0 25]);
title('Put price v.s. Exercise Price');
xlabel('Exercise price');
ylabel('Put price');
legend('EurPut value','內含價值(T=0)');
HW5
function [price, lattice] = LatticeAmCallDivD(S0,X,r,T,sigma,N,D,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
S0a=S0-D*exp(-r*tau*deltaT);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0 , -X+S0a*(u^j)*(d^(N-j)));
end
for i=N-1:-1:tau
for j=0:i
lattice(i+1,j+1) = max( -X+S0a*u^j*d^(i-j) , ...
exp(-r*deltaT) *(p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
end
end
for i=tau-1:-1:0
for j=0:i lattice(i+1,j+1) = max( -X+S0a*u^j*d^(i-j)+D*exp(-r*(tau-i)*deltaT) , ...
exp(-r*deltaT) *(p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
end
end
price = lattice(1,1);
function [price, lattice] = LatticeAmCallDivP(S0,X,r,T,sigma,N,div,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0 , S0*(u^j)*(d^(N-j))*(1-div)-X);
end
for i=N-1:-1:tau
for j=0:i
lattice(i+1,j+1) = max(-X+S0*u^j*d^(i-j)*(1-div) , ...
exp(-r*deltaT) *(p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
end
end
for i=tau-1:-1:0
for j=0:i
lattice(i+1,j+1) = max(-X+S0*u^j*d^(i-j) , ...
exp(-r*deltaT) *(p * lattice(i+2,j+2) + (1-p) * lattice(i+2,j+1)));
end
end
price = lattice(1,1);
function [price, lattice] = LatticeEurCallDivD(S0,X,r,T,sigma,N,D,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
S0a=S0-D*exp(-r*tau*deltaT);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0 , -X+S0a*(u^j)*(d^(N-j)));
end
for i=N-1:-1:0
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
price = lattice(1,1);
function [price, lattice] = LatticeEurCallDivP(S0,X,r,T,sigma,N,div,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0,-X+S0*(u^j)*(d^(N-j))*(1-div));
end
for i=N-1:-1:tau
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*(1-div)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
for i=tau-1:-1:0
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
price = lattice(1,1);
function [price, lattice] = LatticeEurPutDivD(S0,X,r,T,sigma,N,D,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
S0a=S0-D*exp(-r*tau*deltaT);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0 , X-S0a*(u^j)*(d^(N-j)));
end
for i=N-1:-1:0
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
price = lattice(1,1);
function [price, lattice] = LatticeEurPutDivP(S0,X,r,T,sigma,N,div,tau)
deltaT = T/N;
u=exp(sigma * sqrt(deltaT));
d=1/u;
p=(exp(r*deltaT) - d)/(u-d);
lattice = zeros(N+1,N+1);
for j=0:N
lattice(N+1,j+1)=max(0,X-S0*(u^j)*(d^(N-j))*(1-div));
end
for i=N-1:-1:tau
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*(1-div)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
for i=tau-1:-1:0
for j=0:i
lattice(i+1,j+1) = exp(-r*deltaT)*...
(p*lattice(i+2,j+2)+(1-p)*lattice(i+2,j+1));
end
end
price = lattice(1,1);
function [price, lattice] = TriAmPut(S0,X,r,T,sigma,N,lamda)
deltaT = T/N;
u=exp(lamda*sigma * sqrt(deltaT));
d=1/u;
pu=1/(2*lamda^2)+(r-(sigma^2/2))*sqrt(deltaT)/(2*lamda*sigma);
pm=1-1/(lamda^2);
pd=1-pu-pm;
lattice = zeros(N+1,2*N+1);
for j=1:N+1
lattice(N+1,j)=max(0 , -S0*(d^(N-j+1)) + X);
end
for j=N+2:2*N+1
lattice(N+1,j)=max(0 , -S0*(u^(j-N-1)) + X);
end
for i=N-1:-1:0
for j=1:2*i+1
lattice(i+1,j+1) = max( X-S0*u^j*d^(i-j) , ...
(pd * lattice(i+2,j) + pm * lattice(i+2,j+1)+ ...
pu*lattice(i+2,j+2)));
end
end
function [price, lattice] = TriEurPut(S0,X,r,T,sigma,N,lamda)
deltaT = T/N;
u=exp(lamda*sigma * sqrt(deltaT));
d=1/u;
pu=1/(2*lamda^2)+(r-(sigma^2/2))*sqrt(deltaT)/(2*lamda*sigma);
pm=1-1/(lamda^2);
pd=1-pu-pm;
lattice = zeros(N+1,2*N+1);
for j=1:N+1
lattice(N+1,j)=max(0 , -S0*(d^(N-j+1)) + X);
end
for j=N+2:2*N+1
lattice(N+1,j)=max(0 , -S0*(u^(j-N-1)) + X);
end
for i=N-1:-1:0
for j=1:2*i+1
lattice(i+1,j) = exp(-r*deltaT) * ...
(pd * lattice(i+2,j) + pm * lattice(i+2,j+1)+ ...
pu*lattice(i+2,j+2));
end
end
HW6
% BlsMCAVD.m
function [Price, CI] = BlsMCAVD(S0,X,r,T,sigma,NRepl,q)
nuT = ((r-q)-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
Veps = randn(NRepl,1);
Payoff1 = max( 0 , S0*exp(nuT+siT*Veps)-X);
Payoff2 = max( 0 , S0*exp(nuT+siT*(-Veps))-X);
DiscPayoff = exp(-r*T)*0.5*(Payoff1+Payoff2);
[Price, VarPrice, CI] = normfit(DiscPayoff);
% BlsMCAVDP.m
function [Price, CI] = BlsMCAVDP(S0,X,r,T,sigma,NRepl,q)
nuT = ((r-q)-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
Veps = randn(NRepl,1);
Payoff1 = max( 0 ,X-S0*exp(nuT+siT*Veps));
Payoff2 = max( 0 ,X-S0*exp(nuT+siT*(-Veps)));
DiscPayoff = exp(-r*T)*0.5*(Payoff1+Payoff2);
[Price, VarPrice, CI] = normfit(DiscPayoff);
% BlsMCAVP.m
function [Price, CI] = BlsMCAVP(S0,X,r,T,sigma,NRepl)
nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T);
Veps = randn(NRepl,1);
Payoff1 = max( 0 , X-S0*exp(nuT+siT*Veps));
Payoff2 = max( 0 , X-S0*exp(nuT+siT*(-Veps)));
DiscPayoff = exp(-r*T) * 0.5 * (Payoff1+Payoff2);
[Price, VarPrice, CI] = normfit(DiscPayoff);
% BlsMCD.m
function [Price, CI] = BlsMCD(S0,X,r,T,sigma,NRepl,q)
nuT = ((r-q)-0.5*sigma^2)*T;
siT = sigma * sqrt(T);
DiscPayoff = exp(-r*T) * max( 0 , S0*exp(nuT+siT*randn(NRepl,1)) - X);
[Price, VarPrice, CI] = normfit(DiscPayoff);
% BlsMCDP.m
function [Price, CI] = BlsMCDP(S0,X,r,T,sigma,NRepl,q)
nuT = ((r-q)-0.5*sigma^2)*T;
siT = sigma * sqrt(T);
DiscPayoff = exp(-r*T) * max( 0 ,X-S0*exp(nuT+siT*randn(NRepl,1)));
[Price, VarPrice, CI] = normfit(DiscPayoff);
% BlsMCP.m
function [Price, CI] = BlsMCP(S0,X,r,T,sigma,NRepl)
nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T);
DiscPayoff = exp(-r*T) * max( 0 ,X-S0*exp(nuT+siT*randn(NRepl,1)));
[Price, VarPrice, CI] = normfit(DiscPayoff);
% CompBlsMCAVD.m Compare blsprice and BlsMCD200000 and BlsMCAVD100000
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
q=0.05;
NRepl1=100000;
NRepl2=200000;
Bls=blsprice(S0,X,r,T,sigma,q);
randn('seed',0);
[MC200000, CI1] = BlsMCD(S0,X,r,T,sigma,NRepl2,q);
randn('seed',0);
[MCAV100000, CI2] = BlsMCAVD(S0,X,r,T,sigma,NRepl1,q);
Bls
MC200000
CI1
MCAV100000
CI2
% CompBlsMCAVD.m Compare blsprice and BlsMCD200000 and BlsMCAVD100000
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
q=0.05;
NRepl1=100000;
NRepl2=200000;
[BlsC, BlsP]=blsprice(S0,X,r,T,sigma,q);
randn('seed',0);
[MC200000, CI1] = BlsMCDP(S0,X,r,T,sigma,NRepl2,q);
randn('seed',0);
[MCAV100000, CI2] = BlsMCAVDP(S0,X,r,T,sigma,NRepl1,q);
BlsP
MC200000
CI1
MCAV100000
CI2
% ComBlsMCAVP.m Compare blsprice and BlsMc200000 and BlsMCAV100000
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NRepl1=100000;
NRepl2=200000;
[BlsC, BlsP]=blsprice(S0,X,r,T,sigma);
randn('seed',0);
[MCP200000, CI1] = BlsMCP(S0,X,r,T,sigma,NRepl2);
randn('seed',0);
[MCAVP100000, CI2] = BlsMCAVP(S0,X,r,T,sigma,NRepl1);
BlsP
MCP200000
CI1
MCAVP100000
CI2
% ComBlsMCP.m Compare blsprice and BlsMCP1000 and BlsMCP200000
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NRepl1=1000;
NRepl2=200000;
[BlsC, BlsP]=blsprice(S0,X,r,T,sigma);
randn('seed',0);
[MCP1000, CI1000] = BlsMCP(S0,X,r,T,sigma,NRepl1);
randn('seed',0);
[MCP200000, CI200000] = BlsMCP(S0,X,r,T,sigma,NRepl2);
BlsP
MCP1000
CI1000
MCP200000
CI200000
HW8
%BlsHaltonEurCall.m
function Price = BlsHaltonEurCall(S0,X,r,T,sigma,NPoints,Base1,Base2)
nuT = (r-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
rand('seed',0);
H1 = GetHalton(ceil(NPoints/2),Base1);
H2 = GetHalton(ceil(NPoints/2),Base2);
VLog = sqrt(-2*log(H1));
Norm1 = VLog.*cos(2*pi*H2);
Norm2 = VLog.*sin(2*pi*H2);
Norm = [Norm1;Norm2];
DiscPayoff = exp(-r*T) * max( 0 , S0*exp(nuT+siT*Norm) -X );
Price = mean(DiscPayoff);
%{
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NPoints=5000;
%}
%BlsHaltonEurPut.m
function Price = BlsHaltonEurPut(S0,X,r,T,sigma,NPoints,Base1,Base2)
nuT = (r-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
rand('seed',0);
H1 = GetHalton(ceil(NPoints/2),Base1);
H2 = GetHalton(ceil(NPoints/2),Base2);
VLog = sqrt(-2*log(H1));
Norm1 = VLog.*cos(2*pi*H2);
Norm2 = VLog.*sin(2*pi*H2);
Norm = [Norm1;Norm2];
DiscPayoff = exp(-r*T) * max( 0 , X-S0*exp(nuT+siT*Norm));
Price = mean(DiscPayoff);
%{
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NPoints=5000;
%}
function Price = BlsRandBox(S0,X,r,T,sigma,NPoints)
nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T);
% Use Box to generate standard normals
rand('seed',0);
H1 = rand(NPoints,1);
H2 = rand(NPoints,1);
VLog = sqrt(-2*log(H1));
Norm1 = VLog.*cos(2*pi*H2);
Norm2 = VLog.*sin(2*pi*H2);
Norm = [Norm1;Norm2];
DiscPayoff = exp(-r*T) * max( 0 , X - S0*exp(nuT+siT*Norm) );
Price = mean(DiscPayoff);
end
%{
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NPoints=5000;
%}
function Price = BlsRandBoxCall(S0,X,r,T,sigma,NPoints)
nuT = (r - 0.5*sigma^2)*T;
siT = sigma * sqrt(T);
% Use Box to generate standard normals
rand('seed',0);
H1 = rand(NPoints,1);
H2 = rand(NPoints,1);
VLog = sqrt(-2*log(H1));
Norm1 = VLog.*cos(2*pi*H2);
Norm2 = VLog.*sin(2*pi*H2);
Norm = [Norm1;Norm2];
DiscPayoff = exp(-r*T) * max( 0 , S0*exp(nuT+siT*Norm) -X );
Price = mean(DiscPayoff);
end
%{
S0=50;
X=52;
r=0.1;
T=5/12;
sigma=0.4;
NPoints=5000;
%}
function Seq = GetHalton(HowMany, Base)
Seq = zeros(HowMany,1);
NumBits = 1+ceil(log(HowMany)/log(Base));
VetBase = Base.^(-(1:NumBits));
WorkVet = zeros(1,NumBits);
for i=1:HowMany
% increment last bit and carry over if necessary
j=1;
ok = 0;
while ok == 0
WorkVet(j) = WorkVet(j)+1;
if WorkVet(j) < Base
ok = 1;
else
WorkVet(j) = 0;
j = j+1;
end
end
Seq(i) = dot(WorkVet,VetBase);
end
function h=HaltonSeq(m,b)
seq=zeros(m,1);
for n=1:m
n0 = n;
h(n) = 0;
f = 1/b;
while (n0 > 0)
n1 = floor(n0/b);
r = n0 - n1*b;
h(n) = h(n)+f*r;
f = f/b;
n0=n1;
end
end
% CompRandHalton
% Rand
for i=1:1000
U1 = rand()
U2 = rand()
VLog = sqrt(-2*log(U1));
X = VLog * cos(2*pi*U2);
Y = VLog * sin(2*pi*U2);
plot(X,Y,'o');
hold on
end
title('rand Box-Muller distribution');
grid on
HW9
function C=compound(s0,K1,K2,r,T1,T2,sigma,q)
sstar=fzero(@(x) blsprice(x,K2,r,T2-T1,sigma,q)-K1,99);
a1=(log(s0/sstar)+(r-q+(sigma^2)/2)*T1)/(sigma*sqrt(T1));
a2=a1-sigma*sqrt(T1);
b1=(log(s0/K2)+(r-q+(sigma^2)/2)*T2)/(sigma*sqrt(T2));
b2=b1-sigma*sqrt(T2);
C=s0*exp(-q*T2)*binorm(a1,b1,T1,T2)-K2*exp(-r*T2)*binorm(a2,b2,T1,T2)-K1*exp(-r*T1)*normcdf(a2);
end
function Q=binorm(a,b,T1,T2)
fun=@(x,y) (1/(2*pi*sqrt(1-T1/T2)))*exp(-(1/(1-T1/T2))*(x.^2-2*sqrt(T1/T2)*x*y+y.^2)/2);
Q=dblquad(fun,-10,a,-10,b);
end
function z=integrnd(x,y)
T1=0.25;
T2=1;
z=(1/(2*pi*sqrt(1-T1/T2)))*exp(-(1/(1-T1/T2))*(x.^2-2*sqrt(T1/T2)*x*y+y.^2)/2);
end
s0=100;
K1=8;
K2=100;
r=0.1;
T1=1;
T2=1.1;
sigma=0.2;
q=0.05;
Cp=[];
Ce1=[];
for T1=1:0.005:(T2)
C=compound(s0,K1,K2,r,T1,T2,sigma,q);
Cp=[Cp;C];
Ce=blsprice(s0,K1+K2,r,T2,sigma,q);
Ce1=[Ce1;Ce];
end
Cp=[Cp;C]
Ce1=[Ce1;Ce]
plot(Cp)
hold on;
plot(Ce1,'r')
xlabel('T1');
ylabel('price');
axis([1 8 1 15]);
HW10
function Price=BlsMCO(v0,u0,r,qv,qu,T,sigmav,sigmau,N,rho)
nuvT=(r-qv-0.5*sigmav^2)*T;
nuuT=(r-qu-0.5*sigmau^2)*T;
sivT=sigmav*sqrt(T);
siuT=sigmau*sqrt(T);
[x1,x2]=cholesky(rho,N);
DiscPayoff=exp(-r*T)*max(0,v0*exp(nuvT+sivT*x1)-u0*exp(nuuT+siuT*x2));
Price=normfit(DiscPayoff);
function [x1,x2]=cholesky(rho,N)
randn('seed',0);
x1=randn(N,1);
x2=rho*x1+sqrt(1-rho^2)*randn(N,1);
end
function C=outperform(v0,u0,qv,qu,sigmav,sigmau,T,rho)
sigma=sqrt(sigmav^2+sigmau^2-2*rho*sigmav*sigmau);
d1=(log(v0/u0)+(qu-qv+(sigma^2)/2)*T)/(sigma*sqrt(T));
d2=d1-sigma*sqrt(T);
C=v0*exp(-qv*T)*normcdf(d1)-u0*exp(-qu*T)*normcdf(d2);
end
v0=30;
u0=50;
qv=0.2;
qu=0.3;
sigmav=0.25;
sigmau=0.2;
T=1;
rho=0.2;
r=0.1;
BlsO=outperform(v0,u0,qv,qu,sigmav,sigmau,T,rho)
MCO100=BlsMCO(v0,u0,r,qv,qu,T,sigmav,sigmau,100,rho)
MCO1000=BlsMCO(v0,u0,r,qv,qu,T,sigmav,sigmau,1000,rho)
MCO10000=BlsMCO(v0,u0,r,qv,qu,T,sigmav,sigmau,10000,rho)
MCO100000=BlsMCO(v0,u0,r,qv,qu,T,sigmav,sigmau,100000,rho)
HW12
% AssetPaths1.m
function SPaths=AssetPaths1(S0,mu,sigma,T,NSteps,NRepl)
dt = T/NSteps;
nudt = (mu-0.5*sigma^2)*dt;
sidt = sigma*sqrt(dt);
Increments = nudt + sidt*randn(NRepl, NSteps);
LogPaths = cumsum([log(S0)*ones(NRepl,1) , Increments] , 2);
SPaths = exp(LogPaths);
%computeLookbackMCHalton
S0=50;
X=50;
r=0.1;
T=1;
sigma=0.4;
NSteps=300;
NRepl=1000;
PMC=LookbackMC(S0,r,T,sigma,NRepl,NSteps);
PHalton = LookbackHalton(S0,r,T,sigma,NRepl,NSteps);
PMC
PHalton
% HaltonPaths.m
function SPaths=HaltonPaths(S0,mu,sigma,T,NSteps,NRepl)
dt = T/NSteps;
nudt = (mu-0.5*sigma^2)*dt;
sidt = sigma*sqrt(dt);
NRepl = 2*ceil(NRepl/2); % make sure it's even
% Use Box Muller to generate standard normals
RandMat = zeros(NRepl, NSteps);
seeds = myprimes(2*NSteps);
Base1 = seeds(1:NSteps);
Base2 = seeds((NSteps+1):(2*NSteps));
for i=1:NSteps
H1 = GetHalton(NRepl/2,Base1(i));
H2 = GetHalton(NRepl/2,Base2(i));
VLog = sqrt(-2*log(H1));
Norm1 = VLog .* cos(2*pi*H2);
Norm2 = VLog .* sin(2*pi*H2);
RandMat(:,i) = [Norm1 ; Norm2];
end
Increments = nudt + sidt*RandMat;
LogPaths = cumsum([log(S0)*ones(NRepl,1) , Increments] , 2);
SPaths = exp(LogPaths);
%LookbackOptionHalton
function P = LookbackHalton(S0,r,T,sigma,NRepl,NSteps)
Payoff = zeros(NRepl,1);
Path = HaltonPaths(S0,r,sigma,T,NSteps,NRepl);
SMax=max(Path,[],2);
Payoff = max (0, SMax-Path(NSteps+1));
[P,aux,CI] = normfit(exp(-r*T)*Payoff);
%LookbackOptionMC
function [P,CI] = LookbackMC(S0,r,T,sigma,NRepl,NSteps)
Payoff = zeros(NRepl,1);
for i=1:NRepl
Path = AssetPaths1(S0,r,sigma,T,NSteps,1);
SMax = max(Path);
Payoff(i) = max (0, SMax-Path(NSteps+1));
end
[P,aux,CI]=normfit(exp(-r*T)*Payoff);
function p = myprimes(N)
found = 0;
trynumber = 2;
p = [];
while (found < N)
if isprime(trynumber)
p = [p , trynumber];
found = found + 1;
end
trynumber = trynumber + 1;
end
HW13
% EuPutExpl1.m
function price = EuPutExpl1(S0,X,r,T,sigma,Smax,dS,dt)
% set up grid and adjust increments if necessary
M = round(Smax/dS);
dS = Smax/M;
N = round(T/dt);
dt = T/N;
matval = zeros(M+1,N+1);
vetS = linspace(0,Smax,M+1)';
vetj = 0:N;
veti = 0:M;
% set up boundary conditions
matval(:,N+1) = max(X-vetS,0);
matval(1,:) = X*exp(-r*dt*(N-vetj));
matval(M+1,:) = 0;
% set up coefficients
a = 0.5*dt*(sigma^2*veti - r).*veti;
b = 1- dt*(sigma^2*veti.^2 + r);
c = 0.5*dt*(sigma^2*veti + r).*veti;
% solve backward in time
for j=N:-1:1
for i=2:M
matval(i,j) = a(i)*matval(i-1,j+1) + b(i)*matval(i,j+1)+ ...
c(i)*matval(i+1,j+1);
end
end
% find closest point to S0 on the grid and return price
% possibly with a linear interpolation
idown = floor(S0/dS);
iup = ceil(S0/dS);
if idown == iup
price = matval(idown+1,1);
else
price = matval(idown+1,1) + ...
(S0 –(idown+1)*dS)*(matval(iup+1,1) - matval(iup,1))/dS;
end
HW14
function [price,callprice,vetS,vetT] = EuCallImpl(S0,X,r,T,sigma,Smax,dS,dt)
M=round(Smax/dS);
dS=Smax/M;
N=round(T/dt);
dt=T/N;
callprice=zeros(M+1,N+1);
vetS=linspace(0,Smax,M+1)';
vetj=0:N;
veti=0:M;
vetT=linspace(0,T,N+1)';
callprice(:,N+1)=max(vetS-X,0);
callprice(1,:)=0;
callprice(M+1,:)=Smax-X*exp(-r*dt*(N-vetj));
a=0.5*(r*dt*veti-sigma^2*dt*(veti.^2));
b=1+sigma^2*dt*(veti.^2)+r*dt;
c=-0.5*(r*dt*veti+sigma^2*dt*(veti.^2));
coeff=diag(a(3:M),-1)+diag(b(2:M))+diag(c(2:M-1),1);
[L,U]=lu(coeff);
aux=zeros(M-1,1);
for j=N:-1:1
aux(M-1)=-c(M)*callprice(M+1,j);
%matval(2:M,j)=U\(L\(matval(2:M,j+1)+aux));
callprice(2:M,j)=coeff\(callprice(2:M,j+1)+aux);
end
idown=floor(S0/dS);
iup=ceil(S0/dS);
if idown==iup
price=callprice(idown+1,1);
else
price=callprice(idown+1,1)+(S0-idown*dS)*(callprice(idown+2,1)-callprice(idown+1,1))/dS;
end
function [C,callprice,vetS,vetT]=EurCallExp(S0,X,r,T,sigma,Smax,dS,dt)
M=round(Smax/dS);
dS=Smax/M;
N=round(T/dt);
dt=T/N;
callprice=zeros(M+1,N+1);
vetS=linspace(0,Smax,M+1)';
vetj=0:N;
veti=0:M;
vetT=linspace(0,T,N+1)';
callprice(:,N+1)=max(vetS-X,0);
callprice(1,:)=0;
callprice(M+1,:)=Smax-X*exp(-r*dt*(N-vetj));
a=0.5*dt*(sigma^2*veti-r).*veti;
b=1-dt*(sigma^2*veti.^2+r);
c=0.5*dt*(sigma^2*veti+r).*veti;
for j=N:-1:1
for i=2:M
callprice(i,j)=a(i)*callprice(i-1,j+1)+b(i)*callprice(i,j+1)+c(i)*callprice(i+1,j+1);
end
end
idown=floor(S0/dS);
iup=ceil(S0/dS);
if isequal(idown,iup)
C=callprice(idown+1,1);
else
C=callprice(idown+1,1)+(S0-(idown+1)*dS)*(callprice(iup+1,1)-callprice(iup,1))/dS;
end
S0=50;
X=50;
r=0.1;
T=5/12;
sigma=0.4;
Smax=100;
dS1=0.5;
dS2=1;
dS3=2;
dT=5/2400;
c=blsprice(S0,X,r,T,sigma);
ImpcalldS1=EuCallImpl(S0,X,r,T,sigma,Smax,dS1,dT);
ImpcalldS2=EuCallImpl(S0,X,r,T,sigma,Smax,dS2,dT);
ImpcalldS3=EuCallImpl(S0,X,r,T,sigma,Smax,dS3,dT);
c
ImpcalldS1
ImpcalldS2
ImpcalldS3
[price,callprice,vetS,vetT]=EuCallImpl(S0,X,r,T,sigma,Smax,dS1,dT);
%callprice
mesh(vetT,vetS,callprice);
ylabel('Stock price'); xlabel('Time'); title('European Call Option, Implicit Method');
S0=50;
X=50;
r=0.1;
T=5/12;
sigma=0.3;
Smax=100;
dS=2;
dT=5/1200;
[C,callprice,vetS,vetT]=EurCallExp(S0,X,r,T,sigma,Smax,dS,dT);
%callprice
C
mesh(vetT,vetS,callprice);
ylabel('Stock price'); xlabel('Time'); title('European Call Option, Explicit Method');
