% GBM_Euro_call_spread.m
% Modifies GBM_Euro_call.m to calculate spread option on a pair of GBM assets.


clear

K=1

Nrealiz=10^5
S1=100 % stock price at time 0
S2=100
sig1=0.2
sig2=0.3 % volatility
T=1 % end-time (in years)
ri=0.06 % risk-free interest rate
div1=0.03 % dividend yield
div2=0.04
rho=0.5
N=200 % number of time steps

dt=T/N;
nu1dt=(ri-div1-0.5*sig1^2)*dt;
nu2dt=(ri-div2-0.5*sig2^2)*dt;
sig1sdt=sig1*sqrt(dt);
sig2sdt=sig2*sqrt(dt);
srho=sqrt(1-rho^2);

x1init=log(S1); % initial value for x=ln(S)
x2init=log(S2);
%figure
%hold on
tic

sum_CT=0;
sum_CT2=0;
for i=1:Nrealiz
    x1=zeros(N,1);
    x2=zeros(N,1);
    t=zeros(N,1);
    x1(1)=x1init; % initial condition
    x2(1)=x2init;
    n1=randn(N,1);
    n2=randn(N,1);
    r1=n1;
    r2=rho*n1+srho*n2;
    for n=1:(N-1)
        x1(n+1)=x1(n)+nu1dt+sig1sdt*r1(n);
        x2(n+1)=x2(n)+nu2dt+sig2sdt*r2(n);
        t(n+1)=t(n)+dt;
    end
    S1v=exp(x1);
    S2v=exp(x2);
%    plot(t,Sv,'k-');
    S1T=S1v(N);
    S2T=S2v(N);
    CT=max(0,S1T-S2T-K);
    sum_CT=sum_CT+CT;
    sum_CT2=sum_CT2+CT^2;
end

call_value=sum_CT/Nrealiz*exp(-ri*T)
SD=sqrt( (sum_CT2-sum_CT*sum_CT/Nrealiz)*exp(-2*ri*T)/(Nrealiz-1))
SE=SD/sqrt(Nrealiz)
toc