Prev: fast image comparision
Next: Load in part of variable from .mat file
From: Szilard Szasz- Toth on 3 Aug 2010 09:57
Hello,

I am basically trying to convert an American Option to a European Option by inverting Barone-Adesi and Whaley as described in the Paper "Determining the implied volatility for American options using the QAM" - Kutner (1998). I successfully implemented the code using the Newton-Raphson Method for 2 equations, but because the methods seems to have some problems with certain strike values, I wanted to use Matlab's fsolve function as an alternative.

When using fsolve, however, I never get a solution and the supplied exact partial derivatives for the Jacobian matrix do not match the finite-differencing derivatives from fsolve. The supplied derivatives should be correct however (double checked with values from paper) and the solution of the Newton-Raphson Method converges quickly to a solution. I am not sure what I am doing wrong or missing, so any help would be greatly appreciated. Below are some code excerpt. If the code is hardly readable, I would of course send you my entire code, both the Newton-Raphson and fsolve versions. Thanks in advance.

Code:

Execute with test- data:
[price implVola] = convertAmericanToEuropean(1, 0.04, 80, 100, 0.08, 0.25);
% convertAmericanToEuropean_fsolve(isCall, optionPrice, S, X, r, T)
[price implVola] = convertAmericanToEuropean_fsolve(1, 0.04, 80, 100, 0.08, 0.25);


=====================================================
function [optionPrice, optionVola] = convertAmericanToEuropean_fsolve(isCall, optionPrice, S, X, r, T)

if (isCall)
S_Star_0 = S+10; %Call
Sigma_0 = blkimpv(S, X, r, T, optionPrice, [], [], {'call'});
else
S_Star_0 = S-10; %Put
Sigma_0 = blkimpv(S, X, r, T, optionPrice, [], [], {'put'});
end

initialGuess = [S_Star_0; Sigma_0^2]; % Make a starting guess at the solution
options = optimset('Display','iter', 'DerivativeCheck', 'on' , 'Jacobian', 'on');
x = fsolve(@impliedVolatility_BAW_fsolve, initialGuess, options, isCall, optionPrice, S, X, r, 0.0, T);

optionVola = sqrt(x(2));
optionPrice = EurpeanOptionOnCommodity(CallOrPut, S, X, r, 0.0, T, optionVola);
end
=====================================================



=====================================================
function [F,J] = impliedVolatility_BAW_fsolve(x, isCall, marketPrice, S, X, r, b, T)

try

% Set initial starting values
S_Star = x(1);
Sigma = sqrt(x(2));

if (isCall)

% 1.) Set up variables needed to calculate f and g
M = (2*r)/(Sigma^2);
N = (2*b)/(Sigma^2);
K = 1 - exp(-r*T);
q2 = 0.5 * ( -(N-1) + sqrt( (N-1)^2 + (4*M)/K ) );

d1_star = (log(S_Star/X) + (b+0.5*Sigma^2)*T) / (Sigma * sqrt(T));
d2_star = d1_star - Sigma * sqrt(T);
N_d1_Star = normcdf(d1_star,0,1);
%N_d2_Star = normcdf(d2_star,0,1);
n_d1_Star = normpdf(d1_star,0,1);
n_d2_Star = normpdf(d2_star,0,1);

d1 = (log(S/X) + (b+0.5*Sigma^2)*T) / (Sigma * sqrt(T));
d2 = d1 - Sigma * sqrt(T);
%N_d1 = normcdf(d1,0,1);
%N_d2 = normcdf(d2,0,1);
n_d2 = normpdf(d2,0,1);

A2 = (S_Star/q2) * ( 1 - exp((b-r)*T) * N_d1_Star );

c_S_Star = EurpeanOptionOnCommodity(1, S_Star, X, r, b, T, Sigma);
c = EurpeanOptionOnCommodity(1, S, X, r, b, T, Sigma);


% 2.) Calculate f(x) = 0 and g(x) = 0
f = c_S_Star + S_Star/q2 * ( 1 - exp((b-r)*T) * N_d1_Star ) - (S_Star - X);
g = c + A2 * ((S/S_Star)^q2) - marketPrice;

F = [f g];


% 3.) Calculate Jacobian Matrix of f and g

% Substitution variables
y1 = 1 - exp((b-r)*T) * N_d1_Star;
y2 = n_d1_Star / (S_Star * Sigma * sqrt(T));
y3 = exp((b-r)*T);
y4 = b/(Sigma^4) - ( ( (N-1)*b + 2*r/K ) / ( Sigma^4 * sqrt((N-1)^2 + 4*M/K) ) );
y5 = X * sqrt(T) * exp(-r*T) * n_d2 / (2*Sigma);
y6 = -0.5 * ( log(S_Star/X) + b*T ) / ( (Sigma^3 * sqrt(T)) / (4*Sigma) );
y7 = X * sqrt(T) * exp(-r*T) * n_d2_Star / (2*Sigma);

% Partial derivatives of 2x2Jacobi matrix
dg_dS_Star = ((S/S_Star)^q2) * ( (1/q2)*y1 - (S_Star/q2)*y2*y3 - q2*A2*(1/S_Star) );
df_dS_Star = y3 * N_d1_Star + (y1/q2) - (S_Star/q2)*(y2*y3) - 1;
dg_dSigmaSqr = y5 + ((S/S_Star)^q2) * ( A2*log(S/S_Star)*y4 - (S_Star/q2)*y3*y6*n_d1_Star - y1*y4*(S_Star/(q2^2)) );
df_dSigmaSqr = y7 - y1*y4*(S_Star/(q2^2)) - (S_Star/q2)*y3*y6*n_d1_Star;

J = [dg_dS_Star dg_dSigmaSqr ; df_dS_Star df_dSigmaSqr];


else
% Similar code for put
[...]
end
catch
F = [NaN NaN];
J = [NaN NaN ; NaN NaN];
end
end
 | 
Pages: 1
Prev: fast image comparision
Next: Load in part of variable from .mat file