% Euler_Maruyama Method

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);
        if min(w) < 10
            w(101) = 0;
        end
        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);  
    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

eulerc1 = 12*((1+erf(((log(12/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  
   (15*exp(-.05*.5)*(((1+erf(((log(12/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)));
eulerc2 = ((10^2)/(12))*((1+erf(((log(((10^2)/(12))/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  
   (15*exp(-.05*.5)*(((1+erf(((log(((10^2)/(12))/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)));
eulerc3 = ((12/10)^(1-(2*(-.05)/(.25^2))));  %Possibly incorrect division by sigma
actualcallvalue = eulerc1-eulerc2*eulerc3  %closed-form solution (fixed) 
actualcallmatrix = actualcallvalue*actualcallmatrix;

absoluteerror = abs(expectedcallmatrix-actualcallmatrix);
mm9matrix = mean(expectedcallmatrix)*ones(1,100);
mm10matrix = mean(absoluteerror)*ones(1,100);

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

figure()
plot(imatrix,absoluteerror,imatrix,mm10matrix,'r')
axis([0 100 0 0.018])
title('Absolute Error vs. Iteration Number (Euler-Maruyama)')
xlabel('Iteration Number')
ylabel('Absolute Error')


% Milstein Method

expectedcallmatrix2 = zeros(1,100);
actualcallmatrix2 = zeros(1,100);
imatrix = zeros(1,100);

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

    for j = 1:10000 %iterate through at least 10,000 repititions
        [t,w]=milstein([0,.5],12,100);
        if min(w) < 10
            w(101) = 0;
        end
        xvalue2(j) = w(101); %collect the w value at t=.5 for each iteration
        if xvalue2(j)>15
            difference2(j) = xvalue2(j)-15;
        else
            difference2(j) = 0;
        end
        j = j+1;
    end
    meanxvalue2 = mean(difference2);  
    expectedcallvalue2(i) = exp(-.05*.5)*meanxvalue2;  %compute the expected value of the call
    expectedcallmatrix2(i) = expectedcallvalue2(i);
    actualcallmatrix2(i) = 1;
    imatrix(i) = i;
    i=i+1;
end

milsteinc1 = 12*((1+erf(((log(12/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  
   (15*exp(-.05*.5)*(((1+erf(((log(12/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)));
milsteinc2 = ((10^2)/(12))*((1+erf(((log(((10^2)/(12))/15)+(.05+.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)-...  
   (15*exp(-.05*.5)*(((1+erf(((log(((10^2)/(12))/15)+(.05-.5*(.25^2))*(.5))/(.25*sqrt(.5)))/sqrt(2)))/2)));
milsteinc3 = ((12/10)^(1-(2*(-.05)/(.25^2))));  %Possibly incorrect division by sigma
actualcallvalue2 = milsteinc1-milsteinc2*milsteinc3  %closed-form solution (fixed) 
actualcallmatrix2 = actualcallvalue2*actualcallmatrix2;

absoluteerror2 = abs(expectedcallmatrix2-actualcallmatrix2);
mm11matrix = mean(expectedcallmatrix2)*ones(1,100);
mm12matrix = mean(absoluteerror2)*ones(1,100);

figure()
plot(imatrix,expectedcallmatrix2,imatrix,actualcallmatrix2,'r',imatrix,mm11matrix,'g')
title('Expected Call Value vs. Iteration Number with the Actual Call Value Superimposed (Milstein)')
xlabel('Iteration Number')
ylabel('Call Value')

figure()
plot(imatrix,absoluteerror2,imatrix,mm12matrix,'r')
title('Absolute Error vs. Iteration Number (Milstein)')
xlabel('Iteration Number')
ylabel('Absolute Error')


errordifference = absoluteerror-absoluteerror2;  
meanerrormatrix = mean(errordifference)*ones(1,100);  

if mean(errordifference) > 0
    'milstein has less error'
else
    'euler-maruyama has less error'
end

figure()
plot(absoluteerror)
hold on
plot(absoluteerror2,'r')
title('Plot of Absolute Errors (Blue = Euler-Maruyama, Red = Milstein)')
ylabel('Error')
xlabel('Iteration Number')
hold off

figure()
plot(imatrix,errordifference,imatrix,meanerrormatrix,'r')
title('Relative Error Between Euler-Maruyama and Milstein Methods')
ylabel('Error Difference')
xlabel('Iteration Number')