我正在尝试在Matlab中编写一个函数,该函数使用带有矢量输入的Black Scholes公式来计算看涨期权价格。我到目前为止:
function [C] = BlackScholesCall(S,K,t,r,sigma)
%This function calculates the call price per Black-Scholes equation
%INPUT S ... stock price at time 0
% K ... strike price
% r ... interest rate
% sigma ... volatility of the stock price measured as annual standard deviation
% t ... duration in years
%OUTPUT C ... call price
%USAGE BlackScholesCall(S,K,t,r,sigma)
for l = 1:length(K)
for z = 1:length(t)
d1 = (log(S/K(l)) + (r + 0.5*sigma^2)*t(z))/(sigma*sqrt(t(z)));
d2 = d1 - sigma*sqrt(t(z));
N1 = 0.5*(1+erf(d1/sqrt(2)));
N2 = 0.5*(1+erf(d2/sqrt(2)));
C(l) = S*N1-K(l)*exp(-r*t(z))*N2;
end
end
end
F.e。调用我的函数的代码将是
S = 20
K = 16:21
t = 1:1:5
r = 0.02
sigma = 0.25
C = BlackScholesCall(S, K, t, r, sigma)
但是当我将其与Matlab中blsprice函数的结果进行比较时,会得到不同的结果。我怀疑执行循环的方式可能有问题吗?
您得到的结果与,
>> blsprice(S,K,r,t(end),sigma)
ans =
7.1509 6.6114 6.1092 5.6427 5.2102 4.8097
这是因为通过使用C(l) = ...
,您将C
的每个元素覆盖numel(t)
次,因此仅存储/返回z
的每个值的最后计算值。
至少需要使用,
%C(l) = S*N1-K(l)*exp(-r*t(z))*N2;
C(z,l) = S*N1-K(l)*exp(-r*t(z))*N2;
但是您还应该预先分配输出矩阵。也就是说,在任一循环之前,您都应添加
C = nan(numel(K),numel(t));
最后,您应该注意,根本不需要使用任何循环,
[Kmat,tmat] = meshgrid(K,t);
d1 = (log(S./Kmat) + (r + 0.5*sigma^2)*tmat)./(sigma*sqrt(tmat));
d2 = d1 - sigma*sqrt(tmat);
N1 = 0.5*(1+erf(d1/sqrt(2)));
N2 = 0.5*(1+erf(d2/sqrt(2)));
C = S*N1-Kmat.*exp(-r*tmat).*N2;
R版本可能是以下版本。
BlackScholesCall <- function(S, K, tt, r, sigma){
f <- function(.K, .tt){
d1 <- (log(S/.K) + (r + 0.5*sigma^2)*.tt)/(sigma*sqrt(.tt))
d2 <- d1 - sigma*sqrt(.tt)
S*pnorm(d1) - .K*exp(-r*.tt)*pnorm(d2)
}
m <- length(K)
n <- length(tt)
o <- outer(K, tt, f)
last <- if(m > n) o[n:m, n] else o[m, m:n]
c(diag(o), last)
}
BlackScholesCall(S, K, tt, r, sigma)
#[1] 4.703480 4.783563 4.914990 5.059922 5.210161 5.210161 4.809748