How can I optimize code for explicit scheme?

Hi guys,

I am actually trying to price a continuous Asian put option using the explicit scheme. My code gives me a very accurate answer but unfortunately it takes too long to execute. I really need some help in reducing the time significantly. I reduced the time from 11 minutes to about 2 mins but need further reduction. Any help will be really appreciated. Thanks. Here is the code.

function price=solve_asian_put_explicit(r,sigma,K,T,N_time,N_space,S_max,N_v,S0);
%**********************************************************

%**********************************************************
delta_t=T/N_time; %step in time variable
delta_s=(S_max)/N_space; %step in space variable
ind=floor(S0/delta_s);
delta_s=S0/ind;
delta_v=(K*T)/N_v;
rho=delta_t/(delta_s^2);
dt_dv=delta_t/delta_v;
dt_ds=delta_t/(2*delta_s);

%**********************************************************
s=(1:(N_space-1))*delta_s;
s_sqr=s.^2;
v=(0:(N_v+1))*delta_v;
%**********************************************************
phi=zeros(N_space-1,N_v+2);
phi_new=zeros(N_space-1,N_v+2);
%initial condition
for i=1:N_space-1
    for j=1:N_v+2
        phi(i,j)=max((K*T)-v(j),0);
    end
end

%**********************************************************
a=(-r*s*(dt_ds)+(rho*0.5*sigma^2*s_sqr));
b=(1-rho*sigma^2*(s_sqr)-delta_t*r-1.5*s*dt_dv);
c=(2*s*dt_dv);
d=(-0.5*s*dt_dv);
e=(dt_ds*r*s+rho*0.5*sigma^2*s_sqr);
for k=1:N_time
    for i=2:N_space-2
        for j=1:N_v
                phi_new(i,j)= phi(i-1,j)*a(i)+phi(i,j)*b(i)+phi(i,j+1)*c(i)+phi(i,j+2)*d(i)+phi(i+1,j)*e(i);              
        end;
    end;
    phi=phi_new;
end;

%**********************************************************

% return the result
price=phi(ind,1)/T

profile on;
%****************************************
T=1; %maturity
r=0.03; %interest rate
sigma=0.3; %volatility
K=100; %strike price
S0=100; %initial stock price
S_max=S0*exp((r-sigma^2/2)*T+3*sigma*sqrt(T));
N_space=200;
N_time=20000;
N_v=400;

% compute approximate solution
tic

f_approx=solve_asian_put_explicit(r,sigma,K,T,N_time,N_space,S_max,N_v,S0);

time_explicit_scheme=toc
%****************************************

-Cheers

[Walter Roberson]

In your initial condition,

phi(i,j)=max((K*T)-v(j),0);

does not depend upon "i" at all, so why not just do the calculations once for each j and then copy the results by repmat() ?

Answer 1

[Cedric]

You can compute

 A = repmat(a(2:N_space-2), 1, N_v) ;
 B = repmat(b(2:N_space-2), 1, N_v) ;
 C = repmat(c(2:N_space-2), 1, N_v) ;
 D = repmat(d(2:N_space-2), 1, N_v) ;
 E = repmat(e(2:N_space-2), 1, N_v) ;
 for k = 1 : N_time
    phi_new2 = phi(1:N_space-3,1:N_v)   .* A + ...
               phi(2:N_space-2,1:N_v)   .* B + ...
               phi(2:N_space-2,2:N_v+1) .* C + ...
               phi(2:N_space-2,3:N_v+2) .* D + ...
               phi(3:N_space-1,1:N_v)   .* E ;
    phi = phi_new2 ;
 end

which should be more efficient than the double nested loop. Once you are sure that phi_new2 is equivalent to phi_new (check extensively), you make the replacement.

Note that you must use a copy of the original phi when you compute phi_new2 and not the phi that gets out of the triple loop.