r = 0.1; 
a = 0.3;
b = 0.3;
T = 1;
N = 400;      % number of subintervals - needs to be justified!
M = 10000;    % number of simulations  -        "     " 
S0 = 1;
K = 1.1; 

h = T/N;

Z1 = randn(M,N);     % random variables for generating W1
Z2 = randn(M,N);     % random variables for generating W2

c = 0:0.05:1;                % different c's 
tren = zeros(size(c));       % this will hold the answers (for all the c's)
err = zeros(size(c));        % this will hold the stochastic error estimates

for j=1:length(c)

S = S0*ones(M,1);
for i=1:N
    % EM method 
    S = S.*(1 + r*h + a*sqrt(h)*Z1(:,i)  +  b*sqrt(h)*Z2(:,i)./(1+c(j)*S.^2) );
end
ret=exp(-r*T)*(S-K).*(S>K); 
tren(j)=mean(ret);
err(j)=std(ret)/sqrt(M);

end

% plot the results - need to explain them! 
errorbar(c,tren,err)

