exporting patch:
# HG changeset patch
# User parsiad
# Date 1381460384 14400
#      Thu Oct 10 22:59:44 2013 -0400
# Node ID 9e83d2a64ae8a69417fa268feeb672e07828221a
# Parent  7350d2b4eaab53a26cfd933f19a5561a439c1f83
Addition of test cases.

diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blkimpv.m
--- a/inst/blkimpv.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blkimpv.m	Thu Oct 10 22:59:44 2013 -0400
@@ -36,10 +36,23 @@
 ## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
 ## @end deftypefn
 
-function sigma = blkimpv (F, K, r, T, value, limit = 10, tol = 1e-6, option = 1)
+function Volatility = blkimpv (Price, Strike, Rate, Time, Value, ...
+                               Limit = 10, Tolerance = 1e-6, Option = true)
   if (nargin < 5 || nargin > 8)
     print_usage ();
   endif
 
-  sigma = blsimpv (F .* exp (-r .* T), K, r, T, value, limit, 0, tol, option);
+  Volatility = blsimpv (Price .* exp (-Rate .* Time), Strike, Rate, Time, ...
+                        Value, Limit, 0, Tolerance, Option);
 endfunction
+
+## Tests
+%!assert (blkimpv (90, 100, 0.04, 1, [3.4484 13.0563], 1, 1e-6, [1 0]), [0.2 0.2], 1e-4)
+
+## Test input validation
+%!error blkimpv ()
+%!error blkimpv (1)
+%!error blkimpv (1, 2)
+%!error blkimpv (1, 2, 3)
+%!error blkimpv (1, 2, 3, 4)
+%!error blkimpv (1, 2, 3, 4, 5, 6, 7, 8, 2)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blkprice.m
--- a/inst/blkprice.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blkprice.m	Thu Oct 10 22:59:44 2013 -0400
@@ -27,10 +27,26 @@
 ## @seealso{binprice, blsprice}
 ## @end deftypefn
 
-function [C, P] = blkprice (F, K, r, T, sigma)
+function [Call, Put] = blkprice (Price, Strike, Rate, Time, Volatility)
+
   if (nargin != 5)
     print_usage ();
   endif
 
-  [C, P] = blsprice (F .* exp (-r .* T), K, r, T, sigma);
+  [Call, Put] = blsprice (Price .* exp (-Rate .* Time), Strike, Rate, Time, ...
+                          Volatility);
+
 endfunction
+
+## Tests
+%!test
+%! [Call, Put] = blkprice (90:10:110, 100, 0.04, 1, 0.2);
+%! assert (Call, [3.4484 7.6532 13.7316], 1e-4)
+%! assert (Put, [13.0563 7.6532 4.1237], 1e-4)
+
+## Test input validation
+%!error blkprice ()
+%!error blkprice (1)
+%!error blkprice (1, 2)
+%!error blkprice (1, 2, 3)
+%!error blkprice (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsdelta.m
--- a/inst/blsdelta.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsdelta.m	Thu Oct 10 22:59:44 2013 -0400
@@ -31,17 +31,32 @@
 ## @seealso{blsgamma, blslambda, blsprice, blsrho, blstheta, blsvega}
 ## @end deftypefn
 
-function [C, P] = blsdelta (S_0, K, r, T, sigma, q = 0)
+function [CallDelta, PutDelta] = blsdelta (Price, Strike, Rate, Time, ...
+                                           Volatility, Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  d1 = 1 ./ (sigma .* sqrt (T)) .* (log (S_0 ./ K) + (r - q + sigma .^2 / 2) .* T);
+  d1 = 1 ./ (Volatility .* sqrt (Time)) .* (log (Price ./ Strike) + ...
+       (Rate - Yield + Volatility .^2 / 2) .* Time);
 
   phi = normcdf (d1);
 
-  C = exp (-q .* T) .* phi;
-  P = exp (-q .* T) .* (phi - 1);
+  CallDelta = exp (-Yield .* Time) .* phi;
+  PutDelta = exp (-Yield .* Time) .* (phi - 1);
 
 endfunction
+
+## Tests
+%!test
+%! [CallDelta, PutDelta] = blsdelta (90:10:110, 100, 0.04, 1, 0.2, 0.01);
+%! assert (CallDelta, [0.3871 0.5927 0.7586], 1e-4)
+%! assert (PutDelta, [-0.6030 -0.3973 -0.2314], 1e-4)
+
+## Test input validation
+%!error blsdelta ()
+%!error blsdelta (1)
+%!error blsdelta (1, 2)
+%!error blsdelta (1, 2, 3)
+%!error blsdelta (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsgamma.m
--- a/inst/blsgamma.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsgamma.m	Thu Oct 10 22:59:44 2013 -0400
@@ -32,13 +32,25 @@
 ## @seealso{blsdelta, blslambda, blsprice, blsrho, blstheta, blsvega}
 ## @end deftypefn
 
-function G = blsgamma (S_0, K, r, T, sigma, q = 0)
+function Gamma = blsgamma (Price, Strike, Rate, Time, Volatility, Yield = 0)
 
   if (nargin < 4 || nargin > 5)
     print_usage ();
   endif
 
-  sigma_sqrtT = sigma .* sqrt (T);
-  d1 = 1 ./ sigma_sqrtT .* (log (S_0 ./ K) + (r - q + sigma .^2 / 2) .* T);
-  G = exp (-q .* T) .* normpdf (d1) ./ (S_0 .* sigma_sqrtT);
+  sigma_sqrtT = Volatility .* sqrt (Time);
+  d1 = 1 ./ sigma_sqrtT .* (log (Price ./ Strike) + (Rate - Yield + ...
+       Volatility .^2 / 2) .* Time);
+  Gamma = exp (-Yield .* Time) .* normpdf (d1) ./ (Price .* sigma_sqrtT);
+
 endfunction
+
+## Tests
+%assert (blsgamma(90:10:110, 100, 0.04, 1, 0.2, 0.01), [0.0211 0.0191 0.0138], 1e-4)
+
+## Test input validation
+%!error blsgamma ()
+%!error blsgamma (1)
+%!error blsgamma (1, 2)
+%!error blsgamma (1, 2, 3)
+%!error blsgamma (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsimpv.m
--- a/inst/blsimpv.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsimpv.m	Thu Oct 10 22:59:44 2013 -0400
@@ -38,29 +38,36 @@
 ## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
 ## @end deftypefn
 
-function sigma = blsimpv (S_0, K, r, T, value, limit = 10, q = 0, tol = 1e-6, option = true)
+function Volatility = blsimpv (Price, Strike, Rate, Time, Value,        ...
+                               Limit = 10, Yield = 0, Tolerance = 1e-6, ...
+                               Class = true)
+
+  if (nargin < 5 || nargin > 9)
+    print_usage ();
+  endif
 
   ## Get maximum sizes of matrices
-  res_size_row = max ([rows(S_0) rows(K) rows(r) rows(T) rows(value) ...
-                       rows(limit) rows(q) rows(tol) rows(option)]);
-  res_size_col = max ([columns(S_0) columns(K) columns(r) columns(T) ...
-                       columns(value) columns(limit) columns(q) ...
-                       columns(tol) columns(option)]);
+  res_size_row = max ([rows(Price) rows(Strike) rows(Rate) rows(Time)      ...
+                       rows(Value) rows(Limit) rows(Yield) rows(Tolerance) ...
+                       rows(Class)]);
+  res_size_col = max ([columns(Price) columns(Strike) columns(Rate) ...
+                       columns(Time)  columns(Value) columns(Limit) ...
+                       columns(Yield) columns(Tolerance) columns(Class)]);
   res_size = [res_size_row res_size_col];
 
   ## Resize scalars
-  S_0     = resize_scalar (res_size, S_0,    "PRICE");
-  K       = resize_scalar (res_size, K,      "STRIKE");
-  r       = resize_scalar (res_size, r,      "RATE");
-  T       = resize_scalar (res_size, T,      "TIME");
-  value   = resize_scalar (res_size, value,  "VALUE");
-  limit   = resize_scalar (res_size, limit,  "LIMIT");
-  q       = resize_scalar (res_size, q,      "YIELD");
-  tol     = resize_scalar (res_size, tol,    "TOLERANCE");
+  Price     = resize_scalar (res_size, Price,     "PRICE");
+  Strike    = resize_scalar (res_size, Strike,    "STRIKE");
+  Rate      = resize_scalar (res_size, Rate,      "RATE");
+  Time      = resize_scalar (res_size, Time,      "TIME");
+  Value     = resize_scalar (res_size, Value,     "VALUE");
+  Limit     = resize_scalar (res_size, Limit,     "LIMIT");
+  Yield     = resize_scalar (res_size, Yield,     "YIELD");
+  Tolerance = resize_scalar (res_size, Tolerance, "TOLERANCE");
 
-  if (isscalar (option))
-    option = repmat ({option}, res_size)
-  elseif (! isequal (size (option), res_size))
+  if (isscalar (Class))
+    Class = repmat ({Class}, res_size);
+  elseif (! isequal (size (Class), res_size))
     error ("blsimpv: CLASS must be a scalar or cell array of conforming size with other inputs",
            arg_name);
   endif
@@ -68,44 +75,48 @@
   ## Check that all values in the cell array option are either boolean
   ## or the string call, and put. Replace the strings call and put by their
   ## respective boolean and then convert it into a boolean matrix.
-  if (iscell (option))
-    bool_idx = cellfun ("islogical", option);
-    call_idx = strcmpi (option, "call");
-    put_idx  = strcmpi (option, "put");
+  if (iscell (Class))
+    bool_idx = cellfun ("islogical", Class);
+    call_idx = strcmpi (Class, "call");
+    put_idx  = strcmpi (Class, "put");
     if (! all (bool_idx | put_idx | call_idx))
       error ("blsimpv: all values in OPTION must be logical, \"call\", or \"put\"");
     endif
-    option(call_idx) = true;
-    option(put_idx)  = false;
-    option = cell2mat (option);
+    Class(call_idx) = true;
+    Class(put_idx)  = false;
+    Class = cell2mat (Class);
   endif
 
-  sigma = zeros (res_size); # output matrix
+  Volatility = zeros (res_size); # output matrix
 
-  for k = 1:res_size_row
-    for j = 1:res_size_col
+  for j = 1:res_size_row
+    for k = 1:res_size_col
 
       ## Make option set
       ## At the time of writing, fzero did not recognize the
       ## GradObj option.
       opts = optimset ("fzero");
-      opts = optimset (opts, "TolX", tol(j, k), "GradObj", "on");
+      opts = optimset (opts, "TolX", Tolerance(j, k), "GradObj", "on");
 
       ## Construct function
-      if (option(j, k)) # call
-        fun = @(x) blsprice (S_0(j, k), K(j, k), r(j, k), T(j, k), x,
-                             q(j, k)) - value(j, k);
+      if (Class(j, k)) # call
+        fun = @(x) blsprice (Price(j, k), Strike(j, k), Rate(j, k), Time(j, k), x,
+                             Yield(j, k)) - Value(j, k);
       else # put
-        fun = @(x) blsput (S_0(j, k), K(j, k), r(j, k), T(j, k), x,
-                           q(j, k)) - value(j, k);
+        fun = @(x) blsput (Price(j, k), Strike(j, k), Rate(j, k), Time(j, k), x,
+                           Yield(j, k)) - Value(j, k);
       endif
 
       ## Solve
-      [x, fval, exitflag, output] = fzero (fun, [tol limit(j, k)], opts);
-      if (exitflag != 1)
+      try
+        [x, fval, exitflag, output] = fzero (fun, [Tolerance(j, k) Limit(j, k)], opts);
+        if (exitflag != 1)
+          x = NaN;
+        endif
+      catch
         x = NaN;
-      endif
-      sigma(j, k) = x;
+      end_try_catch
+      Volatility(j, k) = x;
 
     endfor
   endfor
@@ -121,12 +132,23 @@
   endif
 endfunction
 
-function [C, v] = blscall (S_0, K, r, T, x, q)
-  [C, P]  = blsprice (S_0, K, r, T, x, q);
-  v       = blsvega  (S_0, K, r, T, x, q);
+function [C, v] = blscall (Price, Strike, Rate, Time, x, Yield)
+  [C, P]  = blsprice (Price, Strike, Rate, Time, x, Yield);
+  v       = blsvega  (Price, Strike, Rate, Time, x, Yield);
 endfunction
 
-function [P, v] = blsput (S_0, K, r, T, x, q)
-  [C, P]  = blsprice (S_0, K, r, T, x, q);
-  v       = blsvega  (S_0, K, r, T, x, q);
+function [P, v] = blsput (Price, Strike, Rate, Time, x, Yield)
+  [C, P]  = blsprice (Price, Strike, Rate, Time, x, Yield);
+  v       = blsvega  (Price, Strike, Rate, Time, x, Yield);
 endfunction
+
+## Tests
+%!assert (blsimpv (100, 100, 0.04, 1, [9.9251 6.0040], 1, 0, 1e-6, [1 0]), [0.2 0.2], 1e-4)
+
+## Test input validation
+%!error blsimpv ()
+%!error blsimpv (1)
+%!error blsimpv (1, 2)
+%!error blsimpv (1, 2, 3)
+%!error blsimpv (1, 2, 3, 4)
+%!error blsimpv (1, 2, 3, 4, 5, 6, 7, 8, 2)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blslambda.m
--- a/inst/blslambda.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blslambda.m	Thu Oct 10 22:59:44 2013 -0400
@@ -14,8 +14,8 @@
 ## along with this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {[@var{CallLambda}, @var{PutLambda}] =} blslambda (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
-## @deftypefnx {Function File} {[@var{CallLambda}, @var{PutLambda}] =} blslambda (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+## @deftypefn  {Function File} {[@var{CallEl}, @var{PutEl}] =} blslambda (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {[@var{CallEl}, @var{PutEl}] =} blslambda (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
 ## Computes elasticity of option under the Black-Scholes model.
 ##
 ## @deftypevar  @var{Price} The current price of the underlying asset.
@@ -33,23 +33,42 @@
 ## There is a bug in MATLAB's version of blslambda in which the deltas of the
 ## option are not discounted by @var{Yield} in the nonzero dividend case. That
 ## is, they compute normcdf(d1) * S / V when they should compute
-## exp(-q*T) * normcdf(d1) * S / V. At the time of writing, this bug is present
-## in the financial toolbox shipped with R2013a. Both this version of blslambda
-## and that shipped with R2013a agree when there are no dividends.
+## exp(-Yield*T) * normcdf(d1) * S / V. At the time of writing, this bug is
+## present in the financial toolbox shipped with R2013a. Both this version of
+## blslambda and that shipped with R2013a agree when there are no dividends.
 ##
 ## @seealso{blsdelta, blsgamma, blsprice, blsrho, blstheta, blsvega}
 ## @end deftypefn
 
-function [C, P] = blslambda (S_0, K, r, T, sigma, q = 0)
+function [CallEl, PutEl] = blslambda (Price, Strike, Rate, Time, ...
+                                      Volatility, Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  [CPrice, PPrice] = blsprice (S_0, K, r, T, sigma, q);
-  [CDelta, PDelta] = blsdelta (S_0, K, r, T, sigma, q);
+  [Call, Put] = blsprice (Price, Strike, Rate, Time, Volatility, Yield);
+  [CallDelta, PutDelta] = blsdelta (Price, Strike, Rate, Time, Volatility, ...
+                                    Yield);
 
-  C = CDelta .* S_0 ./ CPrice;
-  P = PDelta .* S_0 ./ PPrice;
+  CallEl = CallDelta .* Price ./ Call;
+  PutEl  = PutDelta  .* Price ./ Put;
 
 endfunction
+
+## Tests
+%!test
+%! [CallEl, PutEl] = blslambda (90:10:110, 100, 0.04, 1, 0.2);
+%! assert (CallEl, [7.7536 6.2258 5.0647], 1e-4)
+%! assert (PutEl, [-4.8955 -6.3639 -7.8941], 1e-4)
+%!test
+%! [CallEl, PutEl] = blslambda (90:10:110, 100, 0.04, 1, 0.2, 0.01);
+%! assert (CallEl, [7.9108 6.3601 5.1762], 1e-4)
+%! assert (PutEl, [-4.7695 -6.2139 -7.7255], 1e-4)
+
+## Test input validation
+%!error blslambda ()
+%!error blslambda (1)
+%!error blslambda (1, 2)
+%!error blslambda (1, 2, 3)
+%!error blslambda (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsprice.m
--- a/inst/blsprice.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsprice.m	Thu Oct 10 22:59:44 2013 -0400
@@ -32,23 +32,38 @@
 ## @seealso{blskprice, blsdelta, blsgamma, blsimpv, blslambda, blsrho, blstheta, blsvega}
 ## @end deftypefn
 
-function [C, P] = blsprice (S_0, K, r, T, sigma, q = 0)
+function [Call, Put] = blsprice (Price, Strike, Rate, Time, Volatility, ...
+                                 Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  sigma_sqrtT = sigma .* sqrt (T);
+  sigma_sqrtT = Volatility .* sqrt (Time);
 
-  d1 = 1 ./ sigma_sqrtT .* (log (S_0 ./ K) + (r - q + sigma.^2 / 2) .* T);
+  d1 = 1 ./ sigma_sqrtT .* (log (Price ./ Strike) + (Rate - Yield + ...
+       Volatility.^2 / 2) .* Time);
   d2 = d1 - sigma_sqrtT;
 
   phi1 = normcdf (d1);
   phi2 = normcdf (d2);
-  disc = exp (-r .* T);
-  F    = S_0 .* exp ((r - q) .* T);
+  disc = exp (-Rate .* Time);
+  F    = Price .* exp ((Rate - Yield) .* Time);
 
-  C = disc .* (F .*      phi1  - K .*  phi2     );
-  P = disc .* (K .* (1 - phi2) + F .* (phi1 - 1));
+  Call = disc .* (F .*      phi1  - Strike .*  phi2     );
+  Put  = disc .* (Strike .* (1 - phi2) + F .* (phi1 - 1));
 
 endfunction
+
+## Tests
+%!test
+%! [Call, Put] = blsprice (90:10:110, 100, 0.04, 1, 0.2, 0.01);
+%! assert (Call, [4.4037 9.3197 16.1217], 1e-4)
+%! assert (Put, [11.3781 6.3937 3.2952], 1e-4)
+
+## Test input validation
+%!error blsprice ()
+%!error blsprice (1)
+%!error blsprice (1, 2)
+%!error blsprice (1, 2, 3)
+%!error blsprice (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsrho.m
--- a/inst/blsrho.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsrho.m	Thu Oct 10 22:59:44 2013 -0400
@@ -32,18 +32,32 @@
 ## @seealso{blsdelta, blsgamma, blslambda, blsprice, blstheta, blsvega}
 ## @end deftypefn
 
-function [C, P] = blsrho (S_0, K, r, T, sigma, q = 0)
+function [CallRho, PutRho] = blsrho (Price, Strike, Rate, Time, Volatility, Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  d2 = 1 ./ (sigma .* sqrt (T)) .* (log (S_0 ./ K) + (r - q - sigma .^2 / 2) .* T);
+  d2 = 1 ./ (Volatility .* sqrt (Time)) .* (log (Price ./ Strike) + ...
+       (Rate - Yield - Volatility .^2 / 2) .* Time);
 
   phi = normcdf(d2);
+  common = Strike .* Time .* exp (-Rate .* Time);
 
-  common = K .* T .* exp (-r .* T);
-  C = common .* phi;
-  P = common .* (phi - 1);
+  CallRho = common .* phi;
+  PutRho  = common .* (phi - 1);
 
 endfunction
+
+## Tests
+%!test
+%! [CallRho, PutRho] = blsrho (90:10:110, 100, 0.04, 1, 0.2, 0.01);
+%! assert (CallRho, [30.4331 49.9552 67.3271], 1e-4)
+%! assert (PutRho, [-65.6458 -46.1238 -28.7519], 1e-4)
+
+## Test input validation
+%!error blsrho ()
+%!error blsrho (1)
+%!error blsrho (1, 2)
+%!error blsrho (1, 2, 3)
+%!error blsrho (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blstheta.m
--- a/inst/blstheta.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blstheta.m	Thu Oct 10 22:59:44 2013 -0400
@@ -32,27 +32,41 @@
 ## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blsvega}
 ## @end deftypefn
 
-function [C, P] = blstheta (S_0, K, r, T, sigma, q = 0)
+function [CallTheta, PutTheta] = blstheta (Price, Strike, Rate, Time, ...
+                                           Volatility, Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  sqrtT       = sqrt (T);
-  sigma_sqrtT = sigma .* sqrtT;
+  sqrtT       = sqrt (Time);
+  sigma_sqrtT = Volatility .* sqrtT;
 
-  d1 = 1 ./ sigma_sqrtT .* (log (S_0 ./ K) + (r - q + sigma.^2 / 2) .* T);
+  d1 = 1 ./ sigma_sqrtT .* (log (Price ./ Strike) + (Rate - Yield + Volatility.^2 / 2) .* Time);
   d2 = d1 - sigma_sqrtT;
 
   phi1 = normcdf (d1);
   phi2 = normcdf (d2);
 
-  disc  = exp (-q .* T);
-  shift = -disc .* S_0 .* normpdf (d1) .* sigma / 2 ./ sqrtT;
-  t1    = r .* K   .* exp (-r .* T);
-  t2    = q .* S_0 .* disc;
+  disc  = exp (-Yield .* Time);
+  shift = -disc .* Price .* normpdf (d1) .* Volatility / 2 ./ sqrtT;
+  t1    = Rate .* Strike   .* exp (-Rate .* Time);
+  t2    = Yield .* Price .* disc;
 
-  C = shift - t1 .*      phi2  + t2 .*  phi1     ;
-  P = shift + t1 .* (1 - phi2) + t2 .* (phi1 - 1);
+  CallTheta = shift - t1 .*      phi2  + t2 .*  phi1     ;
+  PutTheta  = shift + t1 .* (1 - phi2) + t2 .* (phi1 - 1);
 
 endfunction
+
+## Tests
+%!test
+%! [CallTheta, PutTheta] = blstheta (90:10:110, 100, 0.04, 1, 0.2, 0.01);
+%! assert (CallTheta, [-4.2901 -5.2337 -5.1954], 1e-4)
+%! assert (PutTheta, [-1.3380 -2.3806 -2.4413], 1e-4)
+
+## Test input validation
+%!error blstheta ()
+%!error blstheta (1)
+%!error blstheta (1, 2)
+%!error blstheta (1, 2, 3)
+%!error blstheta (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/blsvega.m
--- a/inst/blsvega.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/blsvega.m	Thu Oct 10 22:59:44 2013 -0400
@@ -32,14 +32,26 @@
 ## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
 ## @end deftypefn
 
-function v = blsvega (S_0, K, r, T, sigma, q = 0)
+function Vega = blsvega (Price, Strike, Rate, Time, Volatility, Yield = 0)
 
   if (nargin < 5 || nargin > 6)
     print_usage ();
   endif
 
-  sqrtT = sqrt (T);
-  d1 = 1 ./ (sigma .* sqrtT) .* (log (S_0 ./ K) + (r - q + sigma.^2 / 2) .* T);
-  v = S_0 .* exp (-q .* T) .* normpdf (d1) .* sqrtT;
+  sqrtT = sqrt (Time);
+  d1 = 1 ./ (Volatility .* sqrtT) .* (log (Price ./ Strike) + ...
+       (Rate - Yield + Volatility.^2 / 2) .* Time);
+  Vega = Price .* exp (-Yield .* Time) .* normpdf (d1) .* sqrtT;
 
 endfunction
+
+## Tests
+%!test
+%assert (blsvega(90:10:110, 100, 0.04, 1, 0.2, 0.01), [34.2115   38.2821   33.3682], 1e-4)
+
+## Test input validation
+%!error blsvega ()
+%!error blsvega (1)
+%!error blsvega (1, 2)
+%!error blsvega (1, 2, 3)
+%!error blsvega (1, 2, 3, 4)
diff -r 7350d2b4eaab -r 9e83d2a64ae8 inst/opprofit.m
--- a/inst/opprofit.m	Tue Oct 08 21:23:55 2013 +0100
+++ b/inst/opprofit.m	Thu Oct 10 22:59:44 2013 -0400
@@ -34,7 +34,7 @@
   elseif (! isbool (PosFlag) && any (PosFlag(:) != 0 & PosFlag(:) != 1))
     error ("opprofit: POSFLAG must be logical, 0 (long) or 1 (short)");
   elseif (! isbool (OptType) && any (OptType(:) != 0 & OptType(:) != 1))
-    error ("opprofit: OPTTYPE must be logical, 0 (long) or 1 (short)");
+    error ("opprofit: OPTTYPE must be logical, 0 (call) or 1 (put)");
   end
 
   pos = -2 * PosFlag + 1;
@@ -42,3 +42,18 @@
   Profit = pos .* (max (opt .* (AssetPrice - Strike), 0) - Cost);
 
 endfunction
+
+## Tests
+%!assert (opprofit (100, 90, 4, 0, 0),  6)
+%!assert (opprofit (80:20:120, 90, 4, 0, 0), [-4 6 26])
+%!assert (opprofit (80:20:120, 80:20:120, 0, 0, 0), [0 0 0])
+%!assert (opprofit (100, 90, 4, [0 0 1 1], [0 1 0 1]), [6 -4 -6 4])
+
+## Test input validation
+%!error opprofit ()
+%!error opprofit (1)
+%!error opprofit (1, 2)
+%!error opprofit (1, 2, 3)
+%!error opprofit (1, 2, 3, 4)
+%!error opprofit (1, 2, 3, 4, 0)
+%!error opprofit (1, 2, 3, 0, 5)
