function [re,op1,opbs]=diffoptionprice(n,So,strike,r,sigma,T)
%[re,op1,opbs]=diffoptionprice(10000,10,8:.1:12,.05,.4,.75)
%estimates the relative error in the BS option price and price under 
% returns distribution which has inverse c.d.f. Finverse. to change from
% logistic change lines. Runs n simulations. the strike price may be a vector.
u=rand(1,n);
%x=log(u./(1-u));   % generates standard logistic
%x=(u<.5).*log(2*u)-(u>.5).*log(2*(1-u)); % generates double exponential
x=log(gaminv(u,2,1));
z=sigma*sqrt(T)*norminv(u)-sigma^2*T/2; % normal returns
c=sigma*sqrt(T/var(x));
mc=mean(exp(c*x));
re=[]; op1=[]; opbs=[];
for i=1:length(strike)
    op1=[op1 mean(max(exp(c*x)*So/mc-exp(-r*T)*strike(i),0))];   % price under F
    opbs=[opbs mean(max(So*exp(z)-exp(-r*T)*strike(i),0))];     % price under BS
end
dif=op1-opbs;
    re=[re dif./(dif+BLSPRICE(So,strike,r,T,sigma,0))];
plot(strike/So,re)
xlabel('Strike price/initial price')
ylabel('relative error in Black Scholes formula')