【正文】 figureplot(Time,w1,39。b39。,Time,w2,39。k39。,Time,Blsprice,39。r:39。)。legend(39。crudeMC39。,39。antitheticMC39。,39。Blsprice39。,2)。title(39。EuropeanoptionpricesusingdifferentMCmethods39。)。xlabel(39。Time39。)。ylabel(39。Optionprices39。)。%%%%Leastsquaresfitoferror=C*ht^q%%%%A=[ones(L_MM,1),log(MM)39。]。b=[log(Error_inf),log(RMS_relative_error)]。22sol=A\b。sol_error=A*solb。resid=[norm(sol_error(:,1)),norm(sol_error(:,2)),norm(sol_error(:,3)),norm(sol_error(:,4))]。qr_inf_c_a_RMSr_c_a=[sol(2,:)。resid]figureloglog(MM,Error_inf(:,1),39。b*39。,MM,Error_inf(:,2),39。39。,MM,RMS_relative_error(:,1),39。go39。,MM,RMS_relative_error(:,2),39。k+39。,MM,MM.^(1/2),39。r:39。)。%axis([10^410^610^410^6])。legend(39。Error\_inf\_c39。,39。Error\_inf\_a39。,39。RMS\_relative\_error\_c39。,39。RMS\_relative\_error\_a39。,39。reference39。)。title(39。loglogplotoftheerror39。)。xlabel(39。M39。)。ylabel(39。Error39。)。子程序：function[w,time]=w_MC(S0,r,sigma,TIME,Z,call_put,K)tic。S_TIME=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*Z)。ifcall_put==1%callPrice=exp(r*TIME).*max(S_TIMEK,0)。else%call_put==2,putPrice=exp(r*TIME).*max(KS_TIME,0)。endw=mean(Price)39。time=toc。function[w,time]=w_MC_antithetic(S0,r,sigma,TIME,Z,call_put,K)tic。Z=Z(1:end/2,:)。TIME=TIME(1:end/2,:)。S_TIME1=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*Z)。S_TIME2=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*(Z))。ifcall_put==1%callPrice1=exp(r*TIME).*max(S_TIME1K,0)。Price2=exp(r*TIME).*max(S_TIME2K,0)。Price=(Price1+Price2)/2。else%call_put==2,putPrice1=exp(r*TIME).*max(KS_TIME1,0)。Price2=exp(r*TIME).*max(KS_TIME2,0)。23Price=(Price1+Price2)/2。endw=mean(Price)39。time=toc。24附錄B：主程序clear。clc。n=100。%numberofreplicationstoevaluateoneoptionm=1e4。%numberofsimulationtrials,eachonebasedonnreplications.Totalreplications=m*n.Sigma=[,]%annualizedassetpricevolatilityRatio=[,]%Ratio=S0/KL_Sigma=length(Sigma)。L_Ratio=length(Ratio)。T=。%timetoexpirationoftheoptionK=100。%strike/exercisepriceoftheoptionr=。%continuouslypoundedriskfreerateofreturnoverthelifeoftheoptionq=0。%dividendyieldcall_put=1。%1forcall。2forputalpha=%degree_of_confidence=1alpharandn(39。state39。,1)。Time_c=zeros(L_Sigma,L_Ratio)。 Time_a=zeros(L_Sigma,L_Ratio)。Time_blsprice=zeros(L_Sigma,L_Ratio)。SE_c=zeros(L_Sigma,L_Ratio)。SE_a=zeros(L_Sigma,L_Ratio)。PERCENT1=zeros(L_Sigma,L_Ratio)。PERCENT2=zeros(L_Sigma,L_Ratio)。BLSPRICE=zeros(L_Sigma,L_Ratio)。foris=1:L_Sigmasigma=Sigma(is)。%annualizedassetpricevolatilityforir=1:L_RatioS0=K*Ratio(ir)。%current/initialpriceoftheunderlyingasset.fprintf(39。Inthiscycle[sigma,S0]=[%d,%d].\n39。,[sigma,S0])%BlackScholesputandcalloptionpricingtic。[Call,Put]=blsprice(S0,K,r,T,sigma,q)。ifcall_put==1%call25Blsprice=Call。else%call_put==2,putBlsprice=Put。endTime_blsprice(is,ir)=tocBLSPRICE(is,ir)=Blsprice。%MonteCarlosimulationsTIME=T*ones(n,m)。%n*mmatrixZ=randn(n,m)。[w1,Time_c(is,ir),SE_c(is,ir)]=w_crude_SE(S0,r,sigma,TIME,Z,call_put,K,n)。[w2,Time_a(is,ir),SE_a(is,ir)]=w_antithetic_SE(S0,r,sigma,TIME,Z,call_put,K,n)。Time_cTime_aPERCENT1(is,ir)=degree_of_confidence(w1,Blsprice,alpha,SE_c(is,ir),m)。PERCENT2(is,ir)=degree_of_confidence(w2,Blsprice,alpha,SE_a(is,ir),m)。endendfprintf(39。displaytheresults\n39。)formatlongBLSPRICEformatSE_cSE_aTime_blspriceTime_cTime_aPERCENT1PERCENT2figureplot(w1,0*ones(m,1),39。b.39。,w2,1*ones(m,1),39。g.39。,[Blsprice,Blsprice],[0,1],39。r*39。)。26legend(39。crudeMC39。,39。antitheticMC39。,39。Blsprice39。,2)。xmin=min([w1。w2])1。xmax=max([w1。w2])+1。ymin=1。ymax=2。axis([xminxmaxyminymax])。title(39。EuropeanoptionpricesusingdifferentMCmethods39。)。xlabel(39。Optionprices39。)。ylabel(39。y39。)。figure[f1,xi1]=ksdensity(w1)。[f2,xi2]=ksdensity(w2)。plot(xi1,f1,39。b.39。,xi2,f2,39。g.39。)。holdon。fmax=max([f1,f2])。plot(Blsprice*ones(2,1),[0,fmax],39。r39。)。holdoff。legend(39。crudeMC39。,39。antitheticMC39。,39。Blsprice39。,2)。title(39。DensityplotsofdifferentMCmethods39。)。xlabel(39。x39。)。ylabel(39。density39。)。子程序：function[w,time,SE_c]=w_crude_SE(S0,r,sigma,TIME,Z,call_put,K,n)tic。S_TIME=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*Z)。ifcall_put==1%callPrice=exp(r*TIME).*max(S_TIMEK,0)。else%call_put==2,putPrice=exp(r*TIME).*max(KS_TIME,0)。endw=mean(Price)39。%anm*1vectorSE_c=std(Price(:))/sqrt(n)。time=toc。function[w,time,SE_a]=w_antithetic_SE(S0,r,sigma,TIME,Z,call_put,K,n)tic。Z=Z(1:end/2,:)。TIME=TIME(1:end/2,:)。S_TIME1=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*Z)。S_TIME2=S0*exp((rsigma^2/2)*TIME+sigma*sqrt(TIME).*(Z))。ifcall_put==1%callPrice1=exp(r*TIME).*max(S_TIME1K,0)。Price2=exp(r*TIME).*max(S_TIME2K,0)。Price=(Price1+Price2)/2。else%call_put==2,put27Price1=exp(r*TIME).*max(KS_TIME1,0)。Price2=exp(r*TIME).*max(KS_TIME2,0)。Price=(Price1+Price2)/2。endw=mean(Price)39。%anm*1vectorSE_a=std(Price(:))/sqrt(n/2)。time=toc；functionpercent=degree_of_confidence(w,Blsprice,alpha,SE,m)z_alpha_2=norminv(1alpha/2,0,1)。percent=sum(abs(wBlsprice)z_alpha_2*SE)/m28