function [m,stderr,mc,stderrc,eff]=asianimportance(r,S0,sig,T,K,k,n,rr)
% Returns:
% m = the value of option calculated using importance sampling
% stderr = standard error of price obtained using a ddifferent measure
% mc = the value of (call asian) option calculated using the original probability measure
% stderrc = standard error for mc
% eff is the eficiency
%input parameters:
%sig = sigma, k=number of time increments in interval [0,T]
% S0=initial stock price
% T= time to maturity
% K = strike price
% n = number of simulations
% rr = the adjustment to the drift of the asset price dynamics
% r =the interest rate parameter. 
%e.g.[m,stderr,mc,stderrc,eff]=asianimportance(0.03,70,0.5,0.25,100,4,10000,2.6)
%assume S0/K is low so that an increase in the drift of the geometric brownian motion is rr
%r=0.03;
%S0=100;


v1=[]; 
FX=[];
GX=[];
mnp=(r-sig^2/2)*T/k;
mn=(r+rr-sig^2/2)*T/k;
sd=sig*sqrt(T/k); 
Y=normrnd(mn,sd,k,n);

%sc= (T/k)*sum(Y-mn)/(sd^2); 
Ycum = cumsum(Y,1);
FX = normpdf(Y,mnp,sd); %a matrix of norm pdfs for f
GX = normpdf(Y,mn,sd); %a matrix of norm pdfs for g
FXprod=prod(FX,1);
GXprod=prod(GX,1);

S = S0*exp(Ycum); 
v1 = exp(-r*T)*max(mean(S,1)-K,0).*(FXprod./GXprod);

m=mean(v1);
stderr=var(v1)/n;

%[mc,stderrc]=asiancrude(r,S0,sig,T,K,k,n);
v2=[]; 


mnp=(r-sig^2/2)*T/k;
Yc=Y-rr*T/k;


%sc= (T/k)*sum(Y-mn)/(sd^2); 
Yc=cumsum(Yc,1);
S = S0*exp(Yc); 
v2= exp(-r*T)*max(mean(S,1)-K,0);

mc=mean(v2);
stderrc=var(v2)/n;
eff=stderrc/stderr;
stderr=sqrt(stderr);
stderrc=sqrt(stderrc);
% disp(['standard errors' num2str(sqrt(var(v1)/n)) num2str(sqrt(var(v1-v2)/n))]);
end