Control variate with Black and Scholes.

[RandomWalk89]

Hey guys,

I wrote this piece of code in order to price a floating strike lookback call option:

<code> % function [Call]=Looback (S,r, sigma,T, M,n)
% MATLAB code for European Lookback Call: Monte-Carlo
%%%%%%%%%% Option parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S = 100; % Value of the underlying
r = 0.1; % Risk free interest rate
sigma = 0.3; % Volatility
T = 0.5;
K=98;
M=100000;
N=200;
Dt=T/N;
% Time to expiry
%%%%%%%%%% Monte-Carlo Method Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% randn(’state’,0)
% M=100000; % Number of Monte-Carlo trials
% n=200; % Set number of steps to compute at in [0,T]
% Compute the time step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dW=sqrt(Dt)*randn(M,N); % Generate array of brownian movements drawn from N(0,dt^2)
W=cumsum(dW,2); % Sum the array cummulativley to obtain Wiener process
t=0t:T; % Array of equal time steps
W=[zeros(M,1),W]; % Set first values to be zero
tt=repmat(t,M,1); % Create matrix for t vals
% Compute the asset path using the solution of the asset path SDE
asset_path=S*exp((r-0.5*sigma^2)*tt+sigma*W);
% Compute the maximum value the asset reaches over the life of the option:
min_vals=min(asset_path,[],2);
% Evaluate the fixed strike lookback call option in each case:
S_final=asset_path,N);
option_values=max(S_final-min_vals,0);
% Discount under risk-neutral assumption
present_vals=exp(-r*T)*option_values;
call_value=mean(present_vals);
call_std = std(present_vals);
confinter = [call_value-1.96*call_std/sqrt(M), call_value+1.96*call_std/sqrt(M)]

Spath = S*cumprod(exp((r-0.5*sigma^2)*Dt+sigma*sqrt(Dt)*randn(M,N)),2);
Spath_final = Spath,N);
mcEC = exp(-r*T)*max(Spath_final-K,0); % Price of a European call with MC

Price_EC = blsprice(S, K, r, T, sigma); % excat Black-Scholes price of European Call

% Control variate
Payoff_control = call_value - mcEC + Price_EC;
P_control_mean = mean(Payoff_control);
P_control_std = std(Payoff_control);
confmc = [P_control_mean-1.96*P_control_std/sqrt(M), P_control_mean+1.96*P_control_std/sqrt(M)] <code>

It works fine except the control variate. I do not get any improvement when implementing it, the confidence interval is slightly the same, [17.1889 17.3805] vs [17.1045 17.3058] with the variance reduction on my last run.

Could you please advise as I'm going nuts?

Cheers guys!!