disp('Simulating the terminal value of option 1000000 times')
U=unifrnd(0,1,1,100000);  
v=fn(U);
disp(['average discounted payoff ' num2str(mean(v))])
Sf=sqrt(var(v));
disp(['Standard error of estimate ' num2str(Sf/sqrt(length(U)) )])
[CALL,PUT] =BLSPRICE(10,10,.05,.25,.2,0) 
disp('***********************************************')
disp('comparing crude MC and BS')
pause
t1=cputime;

ST=[];
	
for i=1:1000 ;
	
S=simulateBS(1432.25,.03,7/252,.195,6);
	
ST=[ST S(7)];
	
end
	disp('option val: simulating whole path 1000 times of GBM for S&P 500 OPTION')
  disp( mean(exp(-.03*7/252)*max((ST-1425),0)))
  t1=[t1 cputime];
   %%%%%%%%%%%%%%%%%%%%%%%%%%%  more efficient
   
	r=.03 ; sigma= .195 ; T=7/252 ; S=1432.25 ;
ST=S*exp(normrnd((r-.5*sigma^2)*T,sigma*sqrt(T),1,100000));
	disp('option val: simulates terminal value 100000 times  for S&P 500 OPTION')

disp(mean(exp(-.03*7/252)*max((ST-1425),0)))
t1=[t1 cputime];
disp(['times    ' num2str(diff(t1))])

disp('BLACK-SCHOLES PRICE USING FORMULA, S&P500 CALL')
disp( BLSPRICE(1432.25,1425,0.03,7/252,0.195,0))
 disp('  *****************************************')
 disp('Call option, S0=10,K=10,r=.05,T=.25,sigma=.2')

pause

t1=cputime;
	
U=rand(1,200000);  % crude Monte Carlo 
	
F=fn(U); % value of 200000 discounted values of the call options
	
m=mean(F) % average discounted value of call = .4743
	
   se=sqrt(var(F)/length(F)) %\ the standard error of the mean
   t1=[t1 cputime];
   printtime(t1)
   


 disp('  *****************************************')
 disp('Using antithetic')
 pause
   t1=[t1 cputime];
	
U=rand(1,100000);  %  use of antithetic random numbers
	
F=.5*(fn(U)+fn(1-U));
mean(F)	
 se=sqrt(var(F)/length(F)) %\ the standard error of the mean
t1=[t1 cputime];
printtime(t1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('  *****************************************')
 disp('Using control')
 pause
t1=cputime;
[est,var1,var2]=control('fn','GG',.438204,200000);
disp(['variance of crude estimator ' num2str(var1/200000)])  % variance of the crude estimator
disp('control variate')
disp([' estimator ' num2str(est) ' variance ' num2str(var2/200000)])
disp(['relative efficiency ' num2str(var1/var2)])
t1=[t1 cputime];
printtime(t1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('  *****************************************')

disp('stratified')
pause
a=.7;
F=a*fn(unifrnd(0,a,1,100000))+(1-a)*fn(unifrnd(a,1,1,100000));
disp([' estimator ' num2str(mean(F)) ' variance ' num2str(var(F)/length(F))])
disp(['relative efficiency ' num2str(var1/(2*var(F)))])
t1=[t1 cputime];
printtime(t1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
disp('  *****************************************')

 
disp('Optimized stratified sample (intervals [0 .6],[.6 .85],[.85 1])')
pause

[est,v,n]=stratified([0 .6 .85 1],200000);
disp([' estimator ' num2str(est) ' variance ' num2str(v)])
disp(['relative efficiency ' num2str(var1/(200000*v))])
t1=[t1 cputime];
printtime(t1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
disp('  *****************************************')
 
disp('combined stratified, antithetic')
pause
U=unifrnd(0,1,1,50000);
F=.37*.5*(fn(.47+.37*U)+fn(.84-.37*U))+.16*.5*(fn(.84+.16*U)+fn(1-.16*U));
disp([' estimator ' num2str(mean(F)) ' variance ' num2str(var(F)/length(F))])
disp(['relative efficiency ' num2str(var1/(4*var(F)))])
t1=[t1 cputime];
printtime(t1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
disp('  *****************************************')

disp('importance sampling')
pause
[est,v]=importance('fn',unifrnd(0,1,1,200000));
disp([' estimator ' num2str(mean(est)) ' variance ' num2str(v/200000)])
disp(['relative efficiency ' num2str(var1/v)])
t1=[t1 cputime];
printtime(t1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
disp('  *****************************************')

disp('Optimal linear combination n=50,000 (time includes all methods)')
pause
plotoptimal
[o,v,b,t1]=optimal(rand(1,50000));
disp(['time   ' num2str(t1)])