求经济金融中的数值方法 试题答案

R语言
#隐式法
library('Matrix')
source('BSput.R')

K<-50       #strike
r<-0.05     #interest
T<-1.0      #1=1 year
S0<-50      #price now
Smax<-3*S0  #max price
sigma<-0.4  #volatility
ds<-1
dt<-0.01
VS<-seq(0,Smax,ds)
VT<-seq(0,T,dt)
M<-length(VS)
N<-length(VT)

M<-M-1
N<-N-1
j<-1:(M+1)

a<- 0.5*r*j*dt - 0.5*sigma^2*j^2*dt
b<- 1+sigma^2*j^2*dt + r*dt
c<- -0.5*r*j*dt - 0.5*sigma^2*j^2*dt

Mat<-bandSparse(M-1,M-1,c(-1,0,1),list(a[2:M],b,c))

f<-matrix(0,M+1,N+1)
f[,N+1]<-matrix(pmax(K-VS,0),M+1,1)
f[1,]<-K
f[M+1,]<-0

ix<-seq(N,1,-1)

for (i in ix)
{
  g<-f[2:M,i+1,drop=F]
  g[1]<-g[1]-a[1]*K
  xx<-solve(Mat,g)
  f[2:M,i]<-as.vector(xx)
}

#prices from BS formula
z<-BSput(VS,0,T,r,sigma,K)
#prices from FDM
y<-f[,1]
x<-VS
plot(x,y,pch=20)
lines(x,z)

#显式法
library('Matrix')
source('BSput.R')

K<-50
r<-0.05
T<-1.0
S0<-50
Smax<-3*S0
sigma<-0.4
ds<-5
dt<-0.01
VS<-seq(0,Smax,ds)
VT<-seq(0,T,dt)
M<-length(VS)
N<-length(VT)
M<-M-1
N<-N-1
j<-1:(M+1)

a<-0.5*r*j*dt - 0.5*sigma^2*j^2*dt
b<-1+sigma^2*j^2*dt + r*dt
c<--0.5*r*j*dt - 0.5*sigma^2*j^2*dt

as<- -a/(1+r*dt)
bs<- (1-sigma^2*j^2*dt)/(1+r*dt)
cs<- -c/(1+r*dt)

Mats<-bandSparse(M+1,M+1,c(0,1,2),list(as,bs,cs))
Mats<-Mats[1:(M-1),]

f<-matrix(0,M+1,N+1)
f[,N+1]<-matrix(pmax(K-VS,0),M+1,1)
f[1,]<-K
f[M+1,]<-0

ix<-seq(N,1,-1)
for (i in ix)
{
  g<-Mats%*%f[,i+1,drop=F]  
  f[2:M,i]<-as.vector(g)  
}

z<-BSput(VS,0,T,r,sigma,K)
y<-f[,1]
x<-VS
plot(x,y,pch=20)
lines(x,z)



#Crank-Nicolsson法
library('Matrix')
source('BSput.R')

K<-50
r<-0.05
T<-1.0
S0<-50
Smax<-3*S0
sigma<-0.4
ds<-1
dt<-0.01
VS<-seq(0,Smax,ds)
VT<-seq(0,T,dt)
M<-length(VS)
N<-length(VT)
M<-M-1
N<-N-1
j<-1:(M+1)

a<- 0.5*r*j*dt - 0.5*sigma^2*j^2*dt
b<- 2+sigma^2*j^2*dt + 2*r*dt
c<- -0.5*r*j*dt - 0.5*sigma^2*j^2*dt

Mat<-bandSparse(M-1,M-1,c(-1,0,1),list(a[2:M],b,c))

as<- -a
bs<- 2-sigma^2*j^2*dt
cs<- -c

Mats<-bandSparse(M-1,M-1,c(-1,0,1),list(as[2:M],bs,cs))

f<-matrix(0,M+1,N+1)
f[,N+1]<-matrix(pmax(K-VS,0),M+1,1)
f[1,]<-K
f[M+1,]<-0

ix<-seq(N,1,-1)
for (i in ix)
{
  g<-Mats%*%f[2:M,i+1,drop=F]
  g[1]<-g[1]+as[1]*K-a[1]*K

  x<-solve(Mat,g)
  f[2:M,i]<-as.vector(x)
}

z<-BSput(VS,0,T,r,sigma,K)
y<-f[,1]
x<-VS
plot(x,y,pch=20)
lines(x,z)

BSput <- function(x,t,T,r,sigma,K)
{
  d2 <- (log(x/K) + (r - 0.5 * sigma^2) * (T - t))/(sigma * sqrt(T - t))
  d1 <- d2 + sigma * sqrt(T - t)
  K * exp(-r * (T - t)) * pnorm(-d2) - x * pnorm(-d1)
}


/******************华丽的分割线******************/
Matlab语言
%隐式法
K=50;
r=0.05;
T=1.0;
S0=50;
Smax=3*S0;
sigma=0.4;
ds=1;
dt=0.01;
VS=0:ds:Smax;
VT=0:dt:T;
M=length(VS);
N=length(VT);

M=M-1;
N=N-1;
j=1:M-1;

a=0.5*r*j*dt - 0.5*sigma^2*j.^2*dt;
b=1+sigma^2*j.^2*dt + r*dt;
c=-0.5*r*j*dt - 0.5*sigma^2*j.^2*dt;
ma=diag(a(2:end),-1);
mc=diag(c(1:end-1),1);
mb=diag(b,0);
Mat=ma+mb+mc;

f=zeros(M+1,N+1);
f(:,end)=max(K-VS,0);
f(1,:)=K;
f(M+1,:)=0;

for i=N:-1:1
    g=f(2:M,i+1);
    g(1)=g(1)-a(1)*K;
    g=Mat\g;
    f(2:M,i)=g;
end
plot(VS,f(:,1),'.-')

%显式法
K=50;
r=0.05;
T=1.0;
S0=50;
Smax=3*S0;
sigma=0.4;
ds=5;
dt=0.01;
VS=0:ds:Smax;
VT=0:dt:T;
M=length(VS);
N=length(VT);

M=M-1;
N=N-1;
j=1:M+1;

as=-(0.5*r*j*dt - 0.5*sigma^2*j.^2*dt)/(1+r*dt);
bs=(1-sigma^2*j.^2*dt)/(1 + r*dt);
cs=-(-0.5*r*j*dt - 0.5*sigma^2*j.^2*dt)/(1+r*dt);
mas=diag(as,0);
mcs=diag(cs(1:end-2),2);
mbs=diag(bs(1:end-1),1);
Mats=mas+mbs+mcs;
Mats=Mats(1:(M-1),:);

f=zeros(M+1,N+1);
f(:,end)=max(K-VS,0);
f(1,:)=K;
f(M+1,:)=0;

for i=N:-1:1
    g=Mats*f(:,i+1);
    f(2:M,i)=g;
end
plot(VS,f(:,1),'.-')

%Crank-Nicolson法
K=50;
r=0.05;
T=1.0;
S0=50;
Smax=3*S0;
sigma=0.4;
ds=1;
dt=0.01;
VS=0:ds:Smax;
VT=0:dt:T;
M=length(VS);
N=length(VT);

M=M-1;
N=N-1;
j=1:M-1;

a=0.5*r*j*dt - 0.5*sigma^2*j.^2*dt;
b=2+sigma^2*j.^2*dt + 2*r*dt;
c=-0.5*r*j*dt - 0.5*sigma^2*j.^2*dt;
ma=diag(a(2:end),-1);
mc=diag(c(1:end-1),1);
mb=diag(b,0);
Mat=ma+mb+mc;

as=-(0.5*r*j*dt - 0.5*sigma^2*j.^2*dt);
bs=(2-sigma^2*j.^2*dt);
cs=-(-0.5*r*j*dt - 0.5*sigma^2*j.^2*dt);
mas=diag(as(2:end),-1);
mcs=diag(cs(1:end-1),1);
mbs=diag(bs,0);
Mats=mas+mbs+mcs;

f=zeros(M+1,N+1);
f(:,end)=max(K-VS,0);
f(1,:)=K;
f(M+1,:)=0;

for i=N:-1:1
    g=Mats*f(2:M,i+1);
    g(1)=g(1)+as(1)*K-a(1)*K;
    x=Mat\g;
    f(2:M,i) = x;
end
plot(VS,f(:,1),'.-')

对Put定价，这是因为Put的边界条件便于处理，Call可以用平价关系获得
显式法不稳定，不要调整显式法的ds和dt我认为Hull书中对Crank-Nicolson的描述是错的，依照书中的逻辑无法得到正确的结果

兄弟，你很厉害，这四道题目是期末作业，需要写好以word形式email给老师的，你能不能把这几题的解题过程发给我啊，谢谢你了。

兄弟，这门课我从来没去上过课，基础几乎为零， 我实在看不懂你写的，并且如果考试不及格要补课的，呵呵，能不能帮我到底啊。