exporting patch:
# HG changeset patch
# User parsiad
# Date 1381197070 14400
#      Mon Oct 07 21:51:10 2013 -0400
# Node ID 83d0c46c16a0a8278b626a864734d905bbb35723
# Parent  7b2b0ca7f368408284f7eac0f4e49f347c382023
Implemented blkimpv, blkprice, blsdelta, blsgamma, blsimpv, blslambda, blsprice, blsrho, blstheta, blsvega, and opprofit functions backwards-engineered from MATLAB Financial Toolbox.

diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blkimpv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blkimpv.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,54 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Volatility} =} blkimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value})
+## @deftypefnx {Function File} {@var{Volatility} =} blkimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit})
+## @deftypefnx {Function File} {@var{Volatility} =} blkimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit}, @var{Tolerance})
+## @deftypefnx {Function File} {@var{Volatility} =} blkimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit}, @var{Tolerance}, @var{Class})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset (a futures contract).
+## @deftypevarx @var{Strike} Exercise price of the futures option.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Value} Price of the European option from which the underlying's volatility is derived.
+## @deftypevarx @var{Limit} (Optional, default = 10) Upper bound of the implied volatility.
+## @deftypevarx @var{Tolerance} (Optional, default = 1e-6) Tolerance with which the root-finding method terminates.
+## @deftypevarx @var{Class} (Optional, default = {'call'}) Option class (call or put). To specify a call option, use a value of true or {'call'}. To specify put options, use a value of false or {'put'}.
+## @end deftypevar
+##
+## Computes the implied volatility under the Black-Scholes model from a given
+## market option price.
+##
+## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
+## @end deftypefn
+
+function sigma = blkimpv(F, K, r, T, value, limit, tol, option)
+
+	if exist('limit') == 0
+		limit = 10;
+	end
+
+	if exist('tol') == 0
+		tol = 1e-6;
+	end
+
+	if exist('option') == 0
+		option = 1;
+	end
+
+	sigma = blsimpv(F .* exp(-r .* T), K, r, T, value, limit, 0, tol, option);
+
+end
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blkprice.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blkprice.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,36 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{Call}, @var{Put}] =} blkprice (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset (a futures contract).
+## @deftypevarx @var{Strike} Exercise price of the futures option.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @end deftypevar
+##
+## Computes the European call and put option prices using the Black-76 model.
+##
+## @seealso{binprice, blsprice}
+## @end deftypefn
+
+function [C, P] = blkprice(F, K, r, T, sigma)
+
+	[C, P] = blsprice(F .* exp(-r .* T), K, r, T, sigma);
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsdelta.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsdelta.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,48 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{CallDelta}, @var{PutDelta}] =} blsdelta (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {[@var{CallDelta}, @var{PutDelta}] =} blsdelta (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the Black-Scholes delta, the rate of change of the option value with
+## respect to the value of the underlying asset.
+##
+## @seealso{blsgamma, blslambda, blsprice, blsrho, blstheta, blsvega}
+## @end deftypefn
+
+function [C, P] = blsdelta(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	d1 = 1 ./ (sigma .* sqrt(T)) .* (log(S_0 ./ K) + (r - q + sigma .^2 / 2) .* T);
+
+	phi = normcdf(d1);
+
+	C = exp(-q .* T) .* phi;
+	P = exp(-q .* T) .* (phi - 1);
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsgamma.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsgamma.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,47 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Gamma} =} blsgamma (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {@var{Gamma} =} blsgamma (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the Black-Scholes gamma, the rate of change of the option delta with
+## respect to the value of the underlying asset.
+##
+## @seealso{blsdelta, blslambda, blsprice, blsrho, blstheta, blsvega}
+## @end deftypefn
+
+function [G] = blsgamma(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	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);
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsimpv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsimpv.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,168 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Volatility} =} blsimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value})
+## @deftypefnx {Function File} {@var{Volatility} =} blsimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit})
+## @deftypefnx {Function File} {@var{Volatility} =} blsimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit}, @var{Yield})
+## @deftypefnx {Function File} {@var{Volatility} =} blsimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit}, @var{Yield}, @var{Tolerance})
+## @deftypefnx {Function File} {@var{Volatility} =} blsimpv (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Value}, @var{Limit}, @var{Yield}, @var{Tolerance}, @var{Class})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Value} Price of the European option from which the underlying's volatility is derived.
+## @deftypevarx @var{Limit} (Optional, default = 10) Upper bound of the implied volatility.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @deftypevarx @var{Tolerance} (Optional, default = 1e-6) Tolerance with which the root-finding method terminates.
+## @deftypevarx @var{Class} (Optional, default = {'call'}) Option class (call or put). To specify a call option, use a value of true or {'call'}. To specify put options, use a value of false or {'put'}.
+## @end deftypevar
+##
+## Computes the implied volatility under the Black-Scholes model from a given
+## market option price.
+##
+## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
+## @end deftypefn
+
+function sigma = blsimpv(S_0, K, r, T, value, limit, q, tol, option)
+
+	if exist('limit') == 0
+		limit = 10;
+	end
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	if exist('tol') == 0
+		tol = 1e-6;
+	end
+
+	if exist('option') == 0
+		option = 1;
+	end
+
+	% Get maximum sizes of matrices
+	res_size_row = max( [size(S_0, 1), size(K, 1), size(r, 1),      ...
+		size(T, 1), size(value, 1), size(limit, 1), size(q, 1), ...
+		size(tol, 1), size(option, 1)] );
+	res_size_col = max( [size(S_0, 2), size(K, 2), size(r, 2),      ...
+		size(T, 2), size(value, 2), size(limit, 2), size(q, 2), ...
+		size(tol, 2), size(option, 2)] );
+	res_size = [res_size_row res_size_col];
+
+	% Resize scalars
+	err = 'Inputs must be scalars or conforming matrices.';
+	if isscalar(S_0)
+		S_0 = S_0 * ones(res_size);
+	elseif isequal(size(S_0), res_size) == 0
+		error(err);
+	end
+	if isscalar(K)
+		K = K * ones(res_size);
+	elseif isequal(size(K), res_size) == 0
+		error(err);
+	end
+	if isscalar(r)
+		r = r * ones(res_size);
+	elseif isequal(size(r), res_size) == 0
+		error(err);
+	end
+	if isscalar(T)
+		T = T * ones(res_size);
+	elseif isequal(size(T), res_size) == 0
+		error(err);
+	end
+	if isscalar(value)
+		value = value * ones(res_size);
+	elseif isequal(size(value), res_size) == 0
+		error(err);
+	end
+	if isscalar(limit)
+		limit = limit * ones(res_size);
+	elseif isequal(size(limit), res_size) == 0
+		error(err);
+	end
+	if isscalar(q)
+		q = q * ones(res_size);
+	elseif isequal(size(q), res_size) == 0
+		error(err);
+	end
+	if isscalar(tol)
+		tol = tol * ones(res_size);
+	elseif isequal(size(tol), res_size) == 0
+		error(err);
+	end
+	if isscalar(option)
+		tmp    = option;
+		option = cell(res_size);
+		for j = 1:res_size_row
+			for k = 1:res_size_col
+				option(j, k) = tmp;
+			end
+		end
+	elseif isequal(size(option), res_size) == 0
+		error(err);
+	end
+
+	% Output matrix
+	sigma = zeros(res_size);
+
+	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');
+
+			% Construct function
+			if isequal(option(j, k), 0) || ...
+					isequal(lower(option(j, k)), {'put'})
+				fun = @(x)   blsput(S_0(j, k), K(j, k), ...
+						r(j, k), T(j, k), x,    ...
+						q(j, k)) - value(j, k);
+			else
+				fun = @(x) blsprice(S_0(j, k), K(j, k), ...
+						r(j, k), T(j, k), x,    ...
+						q(j, k)) - value(j, k);
+			end
+
+			% Solve
+			[x, fval, exitflag, output] = fzero(fun, ...
+					[tol limit(j, k)], opts);
+			if exitflag ~= 1
+				x = NaN;
+			end
+			sigma(j, k) = x;
+
+		end
+	end
+
+end
+
+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);
+end
+
+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);
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blslambda.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blslambda.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,55 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## 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})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the elasticity of an option under the Black-Scholes model.
+## Elasticity measures the percent change in the option price per percent change
+## in the underlying asset price.
+##
+## 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.
+##
+## @seealso{blsdelta, blsgamma, blsprice, blsrho, blstheta, blsvega}
+## @end deftypefn
+
+function [C, P] = blslambda(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	[CPrice, PPrice] = blsprice(S_0, K, r, T, sigma, q);
+	[CDelta, PDelta] = blsdelta(S_0, K, r, T, sigma, q);
+
+	C = CDelta .* S_0 ./ CPrice;
+	P = PDelta .* S_0 ./ PPrice;
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsprice.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsprice.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,54 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{Call}, @var{Put}] =} blsprice (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {[@var{Call}, @var{Put}] =} blsprice (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the European call and put option prices using the Black-Scholes
+## model.
+##
+## @seealso{blskprice, blsdelta, blsgamma, blsimpv, blslambda, blsrho, blstheta, blsvega}
+## @end deftypefn
+
+function [C, P] = blsprice(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	sigma_sqrtT = sigma .* sqrt(T);
+
+	d1 = 1 ./ sigma_sqrtT .* (log(S_0 ./ K) + (r - q + sigma.^2 / 2) .* T);
+	d2 = d1 - sigma_sqrtT;
+
+	phi1 = normcdf(d1);
+	phi2 = normcdf(d2);
+	disc = exp(-r .* T);
+	F    = S_0 .* exp((r - q) .* T);
+
+	C = disc .* (F .*      phi1  - K .*  phi2     );
+	P = disc .* (K .* (1 - phi2) + F .* (phi1 - 1));
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsrho.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsrho.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,49 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{CallRho}, @var{PutRho}] =} blsrho (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {[@var{CallRho}, @var{PutRho}] =} blsrho (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the Black-Scholes rho, the rate of change of the option value with
+## respect to the risk-free interest rate.
+##
+## @seealso{blsdelta, blsgamma, blslambda, blsprice, blstheta, blsvega}
+## @end deftypefn
+
+function [C, P] = blsrho(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	d2 = 1 ./ (sigma .* sqrt(T)) .* (log(S_0 ./ K) + (r - q - sigma .^2 / 2) .* T);
+
+	phi = normcdf(d2);
+
+	common = K .* T .* exp(-r .* T);
+	C = common .* phi;
+	P = common .* (phi - 1);
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blstheta.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blstheta.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,58 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {[@var{CallTheta}, @var{PutTheta}] =} blstheta (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {[@var{CallTheta}, @var{PutTheta}] =} blstheta (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the Black-Scholes theta, the rate of change of the option value with
+## respect to the time-to-expiry.
+##
+## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blsvega}
+## @end deftypefn
+
+function [C, P] = blstheta(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	sqrtT       = sqrt(T);
+	sigma_sqrtT = sigma .* sqrtT;
+
+	d1 = 1 ./ sigma_sqrtT .* (log(S_0 ./ K) + (r - q + sigma.^2 / 2) .* T);
+	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;
+
+	C = shift - t1 .*      phi2  + t2 .*  phi1     ;
+	P = shift + t1 .* (1 - phi2) + t2 .* (phi1 - 1);
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/blsvega.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/blsvega.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,47 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Vega} =} blsvega (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility})
+## @deftypefnx {Function File} {@var{Vega} =} blsvega (@var{Price}, @var{Strike}, @var{Rate}, @var{Time}, @var{Volatility}, @var{Yield})
+##
+## @deftypevar  @var{Price} The current price of the underlying asset.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Rate} The risk-free interest rate.
+## @deftypevarx @var{Time} The time-to-expiry.
+## @deftypevarx @var{Volatility} The volatility of the underlying asset.
+## @deftypevarx @var{Yield} (Optional, default = 0) Annualized, continuously compounded rate of dividends of the underlying asset.
+## @end deftypevar
+##
+## Computes the Black-Scholes vega, the rate of change of the option value with
+## respect to the volatility of the underlying asset.
+##
+## @seealso{blsdelta, blsgamma, blslambda, blsprice, blsrho, blstheta}
+## @end deftypefn
+
+function v = blsvega(S_0, K, r, T, sigma, q)
+
+	if exist('q') == 0
+		q = 0;
+	end
+
+	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;
+
+end
+
diff -r 7b2b0ca7f368 -r 83d0c46c16a0 inst/opprofit.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/inst/opprofit.m	Mon Oct 07 21:51:10 2013 -0400
@@ -0,0 +1,47 @@
+## Copyright (C) 2013 Parsiad Azimzadeh <parsiad.azimzadeh [at] gmail [dot] com>
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU Lesser General Public License as published by the Free
+## Software Foundation; either version 3 of the License, or (at your option) any
+## later version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License
+## for more details.
+##
+## You should have received a copy of the GNU Lesser General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn  {Function File} {@var{Profit} =} opprofit (@var{AssetPrice}, @var{Strike}, @var{Cost}, @var{PosFlag}, @var{OptType})
+##
+## @deftypevar  @var{AssetPrice} The price of the underlying asset at the expiry time.
+## @deftypevarx @var{Strike} The strike price the option is written on.
+## @deftypevarx @var{Cost} The premium paid/charged for the option.
+## @deftypevarx @var{PosFlag} Option position. 0 = long, 1 = short.
+## @deftypevarx @var{OptType} Option type. 0 = call option, 1 = put option.
+## @end deftypevar
+##
+## Computes the profit of an option.
+##
+## @seealso{binprice, blsprice}
+## @end deftypefn
+
+function Profit = opprofit(AssetPrice, Strike, Cost, PosFlag, OptType)
+
+	if any(any(PosFlag ~= 0 & PosFlag ~= 1))
+		error('PosFlag must be 0 (long) or 1 (short)');
+	end
+
+	if any(any(OptType ~= 0 & OptType ~= 1))
+		error('OptType must be 0 (call option) or 1 (put option)');
+	end
+
+	pos = -2 * PosFlag + 1;
+	opt = -2 * OptType + 1;
+
+	Profit = pos .* (max(opt .* (AssetPrice - Strike), 0) - Cost);
+
+end
+
