clear
expectedcallmatrix = zeros(1,100);
actualcallmatrix = zeros(1,100);
imatrix = zeros(1,100);


for i = 1:100  %use a nested loop to obtain different estimated call values
    xvalue = zeros(1,10000); %initialize matrix for w values obtained for each iteration at t=.5
    difference = zeros(1,10000);  %initialize matrix for difference
    expectedcall = zeros(1,10000); %initialize matrix for expected value

    for j = 1:10000 %iterate through at least 10,000 repititions
        [t,w]=euler_maruyama([0,.5],12,100);
        xvalue(j) = w(101); %collect the w value at t=.5 for each iteration
        if xvalue(j)>15
            difference(j) = xvalue(j)-15;
        else
            difference(j) = 0;
        end
        j = j+1;
    end
    meanxvalue = mean(difference); 
    mm(i) = meanxvalue;
    expectedcallvalue(i) = exp(-.05*.5)*meanxvalue;  %compute the expected value of the call
    expectedcallmatrix(i) = expectedcallvalue(i);
    actualcallmatrix(i) = 1;
    imatrix(i) = i;
    i=i+1;
end

actualcallvalue = 12*((1+erf(((log(12/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  %closed-form solution (fixed) 
   (15*exp(-.05*.5)*(((1+erf(((log(12/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)))
actualcallmatrix = actualcallvalue*actualcallmatrix;

mmmatrix = mean(mm)*ones(1,100);

absoluteerror = abs(expectedcallmatrix-actualcallmatrix);
absoluteerroreuler = absoluteerror;  %superfluous

figure()
plot(imatrix,expectedcallmatrix,imatrix,actualcallmatrix,'r',imatrix,mmmatrix,'g')
title('Expected Call Value vs. Iteration Number with the Actual Call Value Superimposed')
xlabel('Iteration Number')
ylabel('Call Value')

figure()
plot(absoluteerror,'r')
title('Absolute Error vs. Iteration Number')
xlabel('Iteration Number')
ylabel('Absolute Error')