-1

I'm trying to write a function in Matlab that calculates the Call price using the Black Scholes formula with vector inputs. I have so far:

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. the code to call my function would be

S = 20
K = 16:21
t = 1:1:5
r = 0.02
sigma = 0.25
C = BlackScholesCall(S, K, t, r, sigma)

But when I compare this with the results of the blsprice function in Matlab, I get different results. I suspect there might be something wrong with the way I did the loop?

bad_coder
  • 11,289
  • 20
  • 44
  • 72
  • Why the tag `r`? You explicitly ask for Matlab help. – Rui Barradas Sep 29 '19 at 20:58
  • Because the differences for the languages in this case are tiny, so if anyone with R knowledge sees this im sure he can answer the question as well. – Question Anxiety Sep 29 '19 at 21:22
  • @QuestionAnxiety: You are asking why Matlab produces a different answer from your code. R is a different language, which makes the notion of diagnosing your code in R ridiculous. Similarly, how any package of R calculates does not predicate how any package of Matlab calculates, or otherwise. In all, R tag is irrelevant to the 2 parts of your question. I don't think it is acceptable to obscure a question -- potentially making some answers obsolete for other users -- for an attempt to fish more views to your question. – Argyll Jul 21 '22 at 16:35

2 Answers2

0

You are getting the same results as,

>> blsprice(S,K,r,t(end),sigma)
ans =
    7.1509    6.6114    6.1092    5.6427    5.2102    4.8097

This is because by using C(l) = ... you are overwriting each element of C numel(t) times, and hence only storing/returning the last calculated values for each value of z.

At a minimum you need to use,

%C(l) = S*N1-K(l)*exp(-r*t(z))*N2;
C(z,l) = S*N1-K(l)*exp(-r*t(z))*N2;

But you should also pre-allocate your output matrix. That is, before either of the loops, you should add

C = nan(numel(K),numel(t));

Finally, you should note that you don't need to use any loops at all,

[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;
Phil Goddard
  • 10,571
  • 1
  • 16
  • 28
0

An R version could be the following.

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
Rui Barradas
  • 70,273
  • 8
  • 34
  • 66