csprice.m

function [callCS,putCS] = csprice(S, X, r, T, sig, skew, kurt, q)
%CSPRICE Corrado and Su (1996) put and call option pricing.
%   Compute European put and call option prices using the Corrado and Su
%   (1996) model.
%
%   [Call,Put] = blsprice(Price, Strike, Rate, Time, Volatility) [Call,Put]
%   = blsprice(Price, Strike, Rate, Time, Volatility, Yield)
%
%   Optional Input: Yield
%
%   Inputs:
%   Price   	- Current price of the underlying asset.
%
%   Strike      - Strike (i.e., exercise) price of the option.
%
%   Rate        - Annualized continuously compounded risk-free rate of
%                 return over the life of the option, expressed as a
%                 positive decimal number.
%
%   Time        - Time to expiration of the option, expressed in years.
%
%   Volatility  - Annualized asset price volatility (i.e., annualized
%                 standard deviation of the continuously compounded asset
%                 return), expressed as a positive decimal number.
%
%   Skewness    - Skewness of the underlying asset returns.
%
%   Kurtosis    - Kurtosis of the underlying asset returns.
%
%   Optional Input:
%   Yield       - Annualized continuously compounded yield of the underlying
%                 asset over the life of the option, expressed as a decimal
%                 number. If Yield is empty or missing. the default value
%                 is zero.
%
%                 For example, this could represent the dividend yield
%                 (annual dividend rate expressed as a percentage of the
%                 price of the security) or foreign risk-free interest rate
%                 for options written on stock indices and currencies,
%                 respectively.
%
%   Outputs:
%   Call        - Price (i.e., value) of a European call option.
%
%   Put         - Price (i.e., value) of a European put option.
%
%   Example:
%      Consider European stock options with an exercise price of $95 that
%      expire in 3 months. Assume the underlying stock pays no dividends,
%      is trading at $100, and has a volatility of 50% per annum, and that
%      the risk-free rate is 10% per annum. Also assume that the skewness 
%      of the underlying equals 0.15 and the kurtosis is equal to 3.15. 
%      Using this data,
%
%      [Call, Put] = csprice(100, 95, 0.1, 0.25, 0.5, 0.15, 3.15) Call =
%         13.7634      
%      Put =
%          6.4179     
%
%
%   Notes:
%   (1) Any input argument may be a scalar or vector. If a scalar, then
%       that value is used to price all options. If more than one input is a
%       vector or matrix, then the dimensions of those non-scalar inputs must
%       be the same.
%   (2) Ensure that Rate, Time, Volatility, and Yield are expressed in
%       consistent units of time.
%
% Semin Ibisevic (2013)
% $Date: 04/29/2013 $


%
% References:
%   Hull, J.C., "Options, Futures, and Other Derivatives", Prentice Hall,
%     5th edition, 2003.
%   Luenberger, D.G., "Investment Science", Oxford Press, 1998.
%   Andreou, P., Charalambous, C., and Martzoukos, S. Pricing and trading
%     european options by combining artificial neural networks and
%     parametric models with implied parameters. European Journal of
%     Operational Research, 185(3):14151433, 2008.
%   Corrado, C.J., and Su T. Skewness and kurtosis in S&P 500 Index returns
%     implied by option prices. Financial Research 19:17592, 1996.


% -------------------------------------------------------------------------
% Input argument checking & default assignment.
if nargin < 7
    error('Finance:csprice:InsufficientInputs', ...
        'Specify Price, Strike, Rate, Time, Volatility, Skewness and Kurtosis')
end

if (nargin < 8) || isempty(q)
    q = zeros(size(S));
end

message = blscheck('blsprice', S, X, r, T, sig, q);
error(message);


% Perform scalar expansion & guarantee conforming arrays.
try
    [S, X, r, T, sig, q] = finargsz('scalar', S, X, r, T, sig, q);
catch
    error('Finance:csprice:InconsistentDimensions', ...
        'Inputs must be scalars or conforming matrices.')
end

% -------------------------------------------------------------------------
% Black and Scholes Price
[call, put] = blsprice(S, X, r, T, sig, q);

% Corrado and Su approximation
d1 = ( log(S./X) + (r-sig).*T + ((sig.*sqrt(T)).^2)./2 ) ./ ( sig.*sqrt(T) );
nd1 = (1./sqrt(2*pi)) * exp(-(1^2)/2);
NNd1 = cdf('Normal',d1,0,1);
Q3 = ( 1/(3*2) ) .* ( S.*exp(-sig.*T) ) .* ( sig .* sqrt(T) ) ...
    .* ( ( 2*sig.*sqrt(T) - d1 ).*nd1 + sig.^2.*T.*NNd1);
Q4 = ( 1/(4*3*2) ) .* ( S.*exp(-sig.*T) ) .* ( sig .* sqrt(T) ) ...
    .* ( ( d1.^2 - 1 - 3.*sig.*sqrt(T) .* (d1 - sig.*sqrt(T) )).*nd1 + ...
    sig.^3.*T.^(3/2).*NNd1);

% Call and Put prices
callCS = call + skew.*Q3 + (kurt-3).*Q4;
if (nargout>1)
    putCS = put + skew.*Q3 + (kurt-3).*Q4;
end