Python Black-Scholes Model and the Basics of Option Pricing

Let St be the price of a stock at time t. The Black-Scholes model is based on the Geometric Brownian Motion (GBM) model which implies that the logarithmic ... Skiptocontent Menu ? PartI:Risk-neutralvaluation,MonteCarlointegrationvs.theBlack-Scholesformula YoucanfindthecodeintheGitHubrepositoryforthisarticle. TheBlack-Scholes(BS)pricingmodelisstilladefactostandardmethodofpricingfinancialoptions. Eventhoughtherehasbeenmuchresearchintoimprovedandpossiblymorerealisticmodels,thefactisthattheBSmodelisimplicitlyassumedinthewaymostoptionpricesarequotedinpractice,intermsofthemodelparametercalledimpliedvolatility.Anyoneworkingwithoptionsinthefinanceindustrywillhavetoknowthebasicsofthisvaluationmethod.Inthistutorial,wewillguideyouthroughtheminimalamountoftheoryneededtoimplementthepricingmethodinPythonandthenlookatabasiccalculationexample. Asyougooverthearticle,feelfreetowatchtheassociatedexplainervideo: TableofContents TheBlack-ScholesModelTheRisk-NeutralWorldExactComputationviatheBlack-ScholesFormula TheBlack-ScholesModel LetSt bethepriceofastockattimet.TheBlack-ScholesmodelisbasedontheGeometricBrownianMotion(GBM)modelwhichimpliesthatthelogarithmicreturnofthestockpricecanbedescribedbyanormaldistributionwhosevarianceisproportionaltothetimestep. ItisalsoapartofthedefinitionoftheGBMthattheselog-returnsarestatisticallyindependentwhenmeasuredoverdisjointtimeintervals.Weassumethattimet ismeasuredinyears,whileSt ismeasuredinthecurrencyinwhichthestockpriceisdenominated,suchasUSDorEUR.μ iscalledtheexpectedreturnandσ isthevolatility.Nowlet’sassumethattheinitialpriceS0=1andlet’ssayexpectedreturnis5% andthevolatilityis20%,inotherwords,μ=0.05 andσ=0.2. Thereasonfortheterm intheexpressionforthelocationparameterofthenormaldistributionisthattheconvexityoftheexponentialintroducesabiaswhichhastobecompensatedfor;thisphenomenonisalsoknownasJensen’sinequality. Thedistributionofthefuturestockpriceoneyearfromnow,giventhatthepricerightnowisequalto1,isknownasthelognormaldistributionandcanbeplottedinPythonusingthescipy.statssubmoduleandthematplotlibpackagetocreatetheactualplot: fromscipy.statsimportlognorm importmatplotlib.pyplotasplt importnumpyasnp importmath mu=0.05 sigma=0.2 S0=1 x=np.linspace(0,2,1000) y=lognorm.pdf(x,sigma,0,math.exp(mu)) plt.plot(x,y) plt.show() Wecanseethatthemostlikelyoutcome,akathe“mode”ormaximumofthedensityfunction,islocatedslightlytotheleftof1,correspondingwiththemathematicalfactthatthemodeofthelognormaldistributioncanbeshowninthiscasetoequal: Thismightdisappointsomepeoplegiventhattheexpectedreturnwassupposedtobe5%,howeverthedistributionisskewedandsothemodeisactuallynotequaltothemeanwhichis: Inotherwords,themostlikelyoutcomeisthatwelose1%ontheinvestment,howevertheskewedprobabilityoflargepositivereturnsimpliesthatonaveragewewillstillmakea5% return! Ifwewanttoseethisstockpriceprocessfromadifferentperspectivewecandoapathsimulation,i.e.randomlygeneratehypotheticaltrajectoriesofthestockpriceaccordingtotheprobabilitylawoftheGBM.WecandothisinPythonjustusingthenumpypackage.Intheexamplebelowwehavesimulated50realizationsofthestockpricepathover1year,dividedinto100uniformtimeincrements: importnumpyasnp importmatplotlib.pyplotasplt Nsim=30 t0=0 t1=1 Nt=100 mu=0.05 sigma=0.2 S0=1 t=np.linspace(t0,t1,Nt) dt=(t1-t0)/Nt S=np.zeros([Nsim,Nt]) S[:,0]=S0 forjinrange(0,Nt-1): S[:,j+1]=S[:,j]*np.exp((mu-sigma**2/2)*dt+sigma*np.sqrt(dt)*np.random.normal(0,1,Nsim)) forjinrange(0,Nsim): plt.plot(t,S[j,:]) plt.show() Theresultsareseeninthenextfigure,itmaylookvaguelylikerealstockpricecharts,atleastifweignorethefactthatrealstockpricessometimesmakeverysudden,sharpdownwardorupwardmovements,astheresultofsomecrisisorotherrandomeventorrevelationwhichsuddenlychangesthemarket’sperceptionofthevalueofastock. Nowlet’sconsideraEuropeancalloptiononthestockinitiatedattimet=0,withanotionalamountof100shares,expirydatet=1 andstrikeprice$1.1.Theoptioncontractgivestheholdertherightbutnottheobligationtobuy100sharesofthestockattheexpirydateoneyearfromnow,forthepricepershareof1dollarand10cents.ThefactthattheoptionhastobeexercisedatthespecificdateofexpiryiswhatismeantbyaEuropeanoption,therearealsoAmericanoptionswhichcanbeexercisedanytimeuntilexpiryandthereareevenothertypeswhichhavedifferentconventionswithregardtoexercise. Whatwouldbetherealizedprofit/lossfromthiscontract,asafunctionofS1,thestockpriceatexpiry?IfIcouldbuythe100sharesfor$110,andthemarketpricehasincreasedtoabove$1.1,ofcourse,Icouldimmediatelyturnaroundandselltheshareagainforthehigherprice.Sothepayoffwouldbe100*S1–1.1.Howeverifthestockpricefellordidn’tincreasebymorethan10%,myresellvaluewouldnotexceedwhatIhadpaidfortheshares,soIwouldnothaveanyinterestinexercisingtheoption.Thepayoffwouldthenbezero.Sothepayoffwouldinanyeventbegivenbytherandomvariable: Wecanplotthisfunctiontovisualizetheasymmetricdependenceofpayoffonthefinaloutcomeofthestockprice: importnumpyasnp importmatplotlib.pyplotasplt k=1.1 defpayoff(x): return100*np.maximum(0,x-k) x=np.linspace(0,2,100) y=payoff(x) plt.plot(x,y) plt.xlabel('Stockpriceatexpiry') plt.ylabel('Payoff') plt.show() TheRisk-NeutralWorld TheBlack-Scholesmodelimpliesthatthevalueofanoptionatatimet priortoexpiryshouldbeequaltotheexpectedpresentvalueofitsfuturepayoff,withjustalittletwist:theexpectationisnotcomputedusingtheactualprobabilitydistributionofthestockprice,evenifwedidactuallybelieveinthemodel’sstatisticalassumptionsabouttheworld.Instead,theexpectationshouldbetakenaccordingtoarisk-neutralprobabilitydistribution,whichmeansthattheexpectedreturnμisreplacedbytherisk-freeinterestrater whilevolatilityisunchanged.Therisk-freeinterestrateistherateofreturnthataninvestorcouldexpecttoreceivebylendingmoneywithoutrunningtheriskoftheborrowerdefaulting;usually,short-termgovernmentbondratesareusedasaproxyforrisk-freeinterestratesbuteventhisassumptionmaybedebatablethesedays.Inourimaginaryrisk-neutralworld,log-returnswouldhavethedistributiongivenby Theoptionpriceattime0wouldthenbeobtainedbycomputingtheexpectation whereEQdenotestherisk-neutralexpectation.Nowlet’ssetthisupinPythonandcomputetheprice,wecanuseMonteCarlointegrationtocomputetheexpectation,meaningthatwedrawalargenumberofrandomsamplesfromthisprobabilitydistribution(whichwouldcorrespondtotheendvaluesofsimulationpathslikethoseweshowedearlier), andcomputethemeantogiveourestimateoftheexpectation.Accordingtothelawoflargenumbersthisestimateapproximatesthetrueexpectationtoarbitraryprecisionifwejustmakethesamplesizelargeenough. importnumpyasnp importmath Nsim=10000 amount_underlying=100 strike=1.1 sigma=0.2 mu=0.06 r=0.015 defpayoff(x): returnamount_underlying*np.maximum(0,x-strike) num0=np.random.normal(0,1,Nsim) S0=15 S1=np.exp(r-sigma**2/2+sigma*num0) C0=math.exp(-r)*np.mean(payoff(S1)) print(C0) NowexecutingthecodeinPythonwegetthefollowingresult: D:\Finxter\Tutorials\Black-Scholes-1>pythonriskneutral.py 4.555089461101134 Whatthismeansinpracticaltermsisthatwithasharepriceof$1,animpliedvolatilitylevelof20%,andarisk-freeinterestrateof1.5%,weshouldexpecttopay$4.555today(plussometransactionfee)foranoptiontobuythe100sharesinoneyearat$1.1pershare. ExactComputationviatheBlack-ScholesFormula ThemethodofMonteCarlointegrationtoobtaintherisk-neutralexpectationofthediscountedpayofffunctionisaverygeneralmethodinthesensethatitwouldworkforallEuropeanoptionsnomatterwhatpayofffunctionweassumed,andwecouldevenhaveexperimentedwithotherassumptionsonthestockpriceprocess.However,inthesimplifiedlognormalworldoftheBSmodel,itturnsoutthatthere’sactuallyaclosed-formequationthatdescribesthecall-optionprice,theso-calledBlack-Scholesformula: whereC0 isthepriceoftheoptionatthestartofthecontract,K isthestrikeprice,t istimetoexpiry, N isthecumulativedistributionfunctionofthestandardnormaldistribution,whiled1 andd2 aregivenby Let’salsotranscribethesefunctionsintoPythoncode: importnumpyasnp fromscipy.statsimportnorm amount_underlying=100 strike=1.1 sigma=0.2 mu=0.06 r=0.015 S0=1 t=1 deffun_d1(sigma,k,t,r,x): return(np.log(x/k)+(r+sigma**2/2)*t)/(sigma*np.sqrt(t)) deffun_d2(sigma,k,t,r,x): returnfun_d1(sigma,k,t,r,x)-sigma*np.sqrt(t) defcall_value(amount_underlying,sigma,k,t,r,x): d1=fun_d1(sigma,k,t,r,x) d2=fun_d2(sigma,k,t,r,x) temp=norm.cdf(d1)*x-norm.cdf(d2)*k*np.exp(-r*t) returnamount_underlying*temp C0=call_value(amount_underlying,sigma,strike,t,r,S0) print(C0) NowifwerunthiswithPython,wegetthefollowingresult: D:\Finxter\Tutorials\Black-Scholes-1>pythonbsformula.py 4.775025500484964 ComparedtoourMonteCarloresultabove,weseethatthere’sadifferenceinthethirddecimal,4.775vs4.777.Weused10000samplesforoursimulation,let’srunitagainwith1000timesthesamplesize,changingtheNsimparameterto10,000,000: D:\Finxter\Tutorials\Black-Scholes-1>pythonriskneutral.py 4.774596150369479 Nowwearegettingclosertotheformula-basedcalculation,indicatingthatthetwodifferentmethodsareindeedequivalent;togetexactlyidenticalanswerswewouldhavetogotothelimitofaninfinitelylargesamplesize.