% Copyright (c) 2007 the authors listed at the following URL, and/or
% the authors of referenced articles or incorporated external code:
% http://en.literateprograms.org/Asian_Option_Pricing_(MATLAB)?action=history&offset=20061105065258
% 
% Permission is hereby granted, free of charge, to any person obtaining
% a copy of this software and associated documentation files (the
% "Software"), to deal in the Software without restriction, including
% without limitation the rights to use, copy, modify, merge, publish,
% distribute, sublicense, and/or sell copies of the Software, and to
% permit persons to whom the Software is furnished to do so, subject to
% the following conditions:
% 
% The above copyright notice and this permission notice shall be
% included in all copies or substantial portions of the Software.
% 
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
% EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
% MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
% IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
% CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
% TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
% SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
% 
% Retrieved from: http://en.literateprograms.org/Asian_Option_Pricing_(MATLAB)?oldid=7915

function [price, pde_sol] = asiancontinuous(S0,K,r,vol,T)
% ASIANCONTINUOUS - implementation of Vecer's PDE Method for Asion Option
% Pricing, minor modifications by Charles-Albert Lehalle.
% use:
% [price, pde_sol] = asiancontinuous(100, 110, .04, .12, .5)
% pde_sol.plot()

%Implementation of Vecer's PDE Method
%Vecer, J. (2002): "Unified Asian Pricing",  Risk, Vol. 15, No. 6, 113-116 

N = 200;       %number of subintervals in space 
M = 200;       %number of subintervals in time

%more points -> higher precision, but slower

%Xmesh  x
xmin = -1;
xmax = 1;
x = linspace(xmin, xmax, N+1);

%Tspan
t = linspace(0, T, M+1);


m = 0;
sol = pdepe(m, @pdef, @pdeic, @pdebc, x, t);
pde_sol = struct('x', x, 't', t, 'u', sol(:,:,1), 'plot', @plot_sol);


%Output of the value of the option
X_0 = (1-exp(-r*T))*S0/r/T - exp(-r*T)*K;
x0 = X_0/S0;

uout = pdeval(m,x,sol(M+1,:),x0);
price = uout*S0;


fprintf( '\n\nPrice of Asian Option is %8.6f\n\n', price);

function [c, f, s] = pdef(x, t, u, DuDx)
c = 1;

    f = 0.5*vol^2*( (1-exp(-r*t))/(r*T) - x )^2*DuDx;
    s = vol^2*((1-exp(-r*t))/(r*T) - x)*DuDx;

end


function u0 = pdeic(x)
u0 = max(x, 0);
end


function [pl, ql, pr, qr] = pdebc(xl, ul, xr, ur, t)
pl = ul;
ql = 0;
pr = ur - xr;
qr = 0;
end


function h = plot_sol
figure('Color',[0.9412 0.9412 0.9412 ]);
surf(pde_sol.x, pde_sol.t, pde_sol.u, 'edgecolor', 'none');
axis([min(pde_sol.x) max(pde_sol.x) min(pde_sol.t) max(pde_sol.t) min(min(pde_sol.u)) max(max(pde_sol.u))]);
xlabel('X');ylabel('t');
end


end