Stochastic Process Models in Finance

Dr Michael Kelly

Wolfram Technology Group

The Ito Stochastic Process forms the basis of all derivative pricing in finance. All options are special cases of the diffusion equation which is solved by the Black-Scholes partial differential equation (PDE) and because of the Feynman-Kac formula this can also be represented by a corresponding stochastic differential equation (SDE) which can now be described by the Mathematica function ItoProcess[]. For instance stocks and indices can be modeled as Geometric Brownian Motion, which is not only a stochastic process function in Mathematica but is also another example of an Ito process. In the same way the ItoProcess[] can be used to describe bonds with stochastic functions like the Brownian Bridge process and stochastic interest rate models can be described by the Ornstein-Uhlenbeck and the Cox-Ingersoll-Ross processes, all of which are now in the finance platform.

StratonovichProcess

A Stratonovich process corresponding to the GeometricBrownianMotionProcess:

StratonovichProcess[GeometricBrownianMotionProcess[μ,σ,x0]]

A Stratonovich process corresponding to the BrownianBridgeProcess:

StratonovichProcess[BrownianBridgeProcess[σ,{(Subscript[t, 1]),a},{(Subscript[t, 2]),b}]]

A Stratonovich process corresponding to the OrnsteinUhlenbeckProcess:

StratonovichProcess[OrnsteinUhlenbeckProcess[μ,σ,θ,x0]]

A Stratonovich process corresponding to the CoxIngersollRossProcess:

StratonovichProcess[CoxIngersollRossProcess[μ,σ,θ,x0]]

Ito process

An Ito process corresponding to the WienerProcess:

ItoProcess[WienerProcess[μ,σ]]

An Ito process corresponding to the GeometricBrownianMotionProcess:

ItoProcess[GeometricBrownianMotionProcess[μ,σ,x0]]

An Ito process corresponding to the BrownianBridgeProcess:

ItoProcess[BrownianBridgeProcess[σ,{(Subscript[t, 1]),a},{(Subscript[t, 2]),b}]]

An Ito process corresponding to the OrnsteinUhlenbeckProcess:

ItoProcess[OrnsteinUhlenbeckProcess[μ,σ,θ,x0]]

An Ito process corresponding to the CoxIngersollRossProcess:

ItoProcess[CoxIngersollRossProcess[μ,σ,θ,x0]]

Compute cross-covariance of the Ornstein-Uhlenbeck process x(t) and its underlying Wiener process w(t):

joint𝒫=ItoProcess[\[DifferentialD]x[t]==θ(μ-x[t])\[DifferentialD]t+σ \[DifferentialD]w[t],{x[t],w[t]},{x,x0},t,w\[Distributed]WienerProcess[]]

Covariance[joint𝒫[t],1,2]

Geometric Brownian Motion

The famous Black–Scholes model for daily stock prices is built upon the notion that such prices behave like geometric Brownian motion. It is because this process when viewed at any fixed time t has the properties of the LogNormal distribution that they were able to integrate values to solve for any option or random function of the LogNormal distribution.

A geometric Brownian motion process is a special ItoProcess:

ItoProcess[GeometricBrownianMotionProcess[μ,σ,a]]

As well as StratonovichProcess:

StratonovichProcess[GeometricBrownianMotionProcess[μ,σ,a]]

GeometricBrownianMotionProcess is not weakly stationary:

WeakStationarity[GeometricBrownianMotionProcess[μ,σ,a]]

The single time slice x(t) follows a LogNormalDistribution

SliceDistribution[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])],t]

Multi-time slice follows a LogMultinormalDistribution:

GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])][{1,2}]

Mean and variance functions:

Mean[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])][t]]//Simplify

Variance[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])][t]]//Simplify

Covariance function:

CovarianceFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])],s,t]

Correlation function:

CorrelationFunction[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])],s,t]

First-order probability density function:

PDF[GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])][t],x]

Compute an expectation:

Expectation[x[t]^2,x\[Distributed]GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])]]

Calculate the probability of an event:

Probability[x[t]<1,x\[Distributed]GeometricBrownianMotionProcess[μ,σ,(Subscript[x, 0])]]

Define a transformed GeometricBrownianMotionProcess:

𝒫=TransformedProcess[t b[t],b\[Distributed]GeometricBrownianMotionProcess[0,2,1],t];

Simulate the process:

data=RandomFunction[𝒫,{0,0.3,0.01},3];

ListLinePlot[data,PlotRange->All]

Mean and variance functions:

{Mean[𝒫[t]],Variance[𝒫[t]]}

Covariance function:

CovarianceFunction[𝒫,t,s]

Using GBM to estimate daily financial data

fdata=FinancialData["SP500",DatePlus[-360]];

DateListPlot[fdata,PlotLabel->"S&P 500",Joined->True, Filling->Axis]

fvalues=Transpose[{Range[0,Length[fdata]-1],Transpose[fdata][[2]]}];

eproc=EstimatedProcess[fvalues,GeometricBrownianMotionProcess[μ,σ,x]]

sims=RandomFunction[eproc,{0,Length[fdata]-1,1},10]

Length[sims["Paths"]]

ListLinePlot[Join[{fvalues},sims["Paths"]],PlotStyle->Prepend[Table[ColorData[1,2],{Length[sims["Paths"]]+1}],Directive[Thick,ColorData[1,1]]],PlotLegends->{"data","process"}]

BrownianBridgeProcess

ItoProcess[BrownianBridgeProcess[σ,{(Subscript[t, 1]),a},{(Subscript[t, 2]),b}]]

Simulate a path of a Brownian bridge process pinned at two different values:

data=RandomFunction[BrownianBridgeProcess[{0,-1},{2,1}],{0,2,.01},3];

ListLinePlot[data["Paths"]]

stochastic exponential function

Define the stochastic exponential function:

𝒫=TransformedProcess[(E^(b[t]-(t/2))),b\[Distributed]WienerProcess[],t];

data=RandomFunction[𝒫,{0,10,0.01},2];

ListLinePlot[data,PlotRange->All]

{Mean[𝒫[t]],Variance[𝒫[t]]}

The corresponding differential equation is du[t] == u[t] db[t]:

𝒫2=ItoProcess[\[DifferentialD]u[t]==u[t]\[DifferentialD]b[t],u[t],{u,1},t,b\[Distributed]WienerProcess[]];

data=RandomFunction[𝒫2,{0,10,0.01},2];

ListLinePlot[data,PlotRange->All]

{Mean[𝒫2[t]],Variance[𝒫2[t]]}

Heston Stochastic Volatility Model

Define the Ito process of the correlated two-dimensional Wiener process:

cW[ρ_]:=ItoProcess[{{0,0},IdentityMatrix[2]},{{w1,w2},{0,0}},t,{{1,ρ},{ρ,1}}];

Define the Heston model by SDEs driven by the correlated 2D Wiener process:

hm=ItoProcess[{
\[DifferentialD]s[t]==μ s[t]\[DifferentialD]t+(Sqrt[r[t]])s[t]\[DifferentialD](Subscript[w, s])[t],
\[DifferentialD]r[t]==κ (θ-r[t])\[DifferentialD]t+ξ (Sqrt[r[t]])\[DifferentialD](Subscript[w, ν])[t]},
{s[t],r[t]},{{s,r},{(Subscript[s, 0]),(Subscript[r, 0])}},t,{(Subscript[w, s]),(Subscript[w, ν])}\[Distributed]cW[ρ]];

Simulate the model using the stochastic Runge–Kutta scheme:

td=BlockRandom[SeedRandom[1988];RandomFunction[hm/.{μ->0,κ->2,θ->1,ξ->1/2,ρ->-1/3,(Subscript[s, 0])->25,(Subscript[r, 0])->1.25},{0,1,0.005},6,Method->"StochasticRungeKutta"]];

Visualize paths:

Row[{ListLinePlot[td["PathComponent",1],PlotLabel->"Price of the asset"],ListLinePlot[td["PathComponent",2],PlotLabel->"Volatility of the asset"]},Spacer[20]]

Ito and Stratonovich Solutions of the Linear Growth Model

Define Ito and Stratonovich processes with the same SDE dx(t)==x(t) (r dt+σ dw(t)):

ipr=ItoProcess[\[DifferentialD]x[t]==x[t](r \[DifferentialD]t+σ \[DifferentialD]w[t]),x[t],{x,x0},t,w\[Distributed]WienerProcess[]]

spr=StratonovichProcess[\[DifferentialD]x[t]==x[t](r \[DifferentialD]t+σ \[DifferentialD]w[t]),x[t],{x,x0},t,w\[Distributed]WienerProcess[]]

Find mean and variance functions for the Ito process:

{Mean[ipr[t]],Variance[ipr[t]]}//Simplify

The mean and variance functions for the Stratonovich process are different:

{Mean[spr[t]], Variance[spr[t]]}//Simplify

When 0<r<((σ^2)/2), the Ito solution almost surely converges to zero, i.e. the large t limit of probability that process value  x(t) does not exceed (1/t) equals 1:

Assuming[0<r<(σ^2)/2&&σ>0&&x0>0,
Limit[Probability[\[FormalX][t]≤1/t,\[FormalX]\[Distributed]ipr],t->∞]]

Confirm it using simulations:

ListLogPlot[RandomFunction[ipr/.{r->1,σ->2,x0->1},{0,20.,0.002},6,Method->"KPS"],PlotRange->All,Joined->True]

When 0<r<((σ^2)/2), the Stratonovich solution, however, almost surely diverges, i.e. the large t limit of the probability that process value x(t) exceeds t equals 1:

Assuming[0<r<(σ^2)/2&&σ>0&&x0>0,
Limit[Probability[\[FormalX][t]>t,\[FormalX]\[Distributed]spr],t->∞]]

Confirm it using simulations:

ListLogPlot[RandomFunction[spr/.{r->1,σ->2,x0->1},{0,20.,0.002},6,Method->"KPS"],PlotRange->All,Joined->True]

Zero Coupon Bonds

Simulate Zero Coupon Bond Using Chen's Model

chenmodel[θ_,ζ_,β_,η_,{r0_,α0_,σ0_}]:=ItoProcess[{
\[DifferentialD]r[t]==(θ-α[t])\[DifferentialD]t+(Sqrt[r[t]])σ[t]\[DifferentialD]w[t],
\[DifferentialD]α[t]==(ζ-α[t])\[DifferentialD]t+(Sqrt[α[t]])σ[t]\[DifferentialD]w[t],
\[DifferentialD]σ[t]==(β-σ[t])\[DifferentialD]t+(Sqrt[σ[t]])η \[DifferentialD]w[t]
},r[t],{{r,α,σ},{r0,α0,σ0}},t,w\[Distributed]WienerProcess[]
]

Visualize some trajectories using the stochastic Runge–Kutta method, correcting for negative interest rates that result from the use of the approximation:

ListLinePlot[RandomFunction[chenmodel[1,2,1.2,1,{1,1,1}],{0,1,0.01},6,Method->{"StochasticRungeKutta","ProjectionFunction"->Function[{t,xvec},Abs[xvec]]}]]

Evaluate zero-coupon bond maturing at time 1 using the Monte Carlo method. The present value of the bond is given by the mean of the process x(t)==exp(-(∫_0^t)r(s)ds):

zcb[θ_,ζ_,β_,η_,{r0_,α0_,σ0_}]=ItoProcess[\[DifferentialD]x[t]==r[t]\[DifferentialD]t,Exp[-x[t]],{x,0},t,r\[Distributed]chenmodel[θ,ζ,β,η,{r0,α0,σ0}]];

td=RandomFunction[zcb[1,2,1.2,1,{1,1,1}],{0,1,0.02},250,Method->{"StochasticRungeKutta","ProjectionFunction"->Function[{t,xvec},Abs[xvec]]}];
Plot[Mean[td[t]]//Evaluate,{t,0,1},AxesOrigin->{0,0},PlotLabel->Row[{"P(0,1)≃",Mean[td[1]]}]]

Stochastic Logistic Growth Model

Define the stochastic differential equation describing the stochastic logistic growth model:

slgm[r_,xmax_,x0_,σ_]:=ItoProcess[\[DifferentialD]x[t]== x[t](1-((x[t])/xmax))(r \[DifferentialD]t+σ \[DifferentialD]w[t]),x[t],{x,x0},t,w\[Distributed]WienerProcess[]]

Deterministic solution of σ->0 is well known:

detsol[r_,xmax_,x0_]:=Function[t,(((E^(r t)) x0 )/(((E^(r t))-1) (x0/xmax)+1))]

Simulate SDE using method of Kloeden–Platen–Schurz of strong order 1.5:

stochsol[r_,xmax_,x0_,σ_]:=RandomFunction[slgm[r,xmax,x0,σ],{0,6,0.005},10,Method->"KloedenPlatenSchurz"]

Visualize stochastic and deterministic solutions:

With[{r=1,xmax=4,x0=1,σ=1/2},Show[ListLinePlot[stochsol[r,xmax,x0,σ]],Plot[detsol[r,xmax,x0][t],{t,0,6},PlotStyle->Directive[Orange,Thick]]]]

KMV-Merton Model

We can use Wolfram Finance Platform’s processes to move to the KMV-Merton model, which assumes a stochastic total value of the firm,

assetValue=1000;defaultBoundary=90; stdevReturns=20;

kmvMERTONFirmValue=ItoProcess[\[DifferentialD]v[t]==μ\[DifferentialD]t+σ \[DifferentialD]w[t],v[t],{v,(Subscript[v, 0])},t,w\[Distributed]WienerProcess[]];

kmvMERTONFirmValueSIMS=RandomFunction[kmvMERTONFirmValue/.{(Subscript[v, 0])->assetValue, μ->defaultBoundary,σ->stdevReturns},{0.,5.,0.001},10];
ListLinePlot[kmvMERTONFirmValueSIMS]

Barath 2004 (Forecasting Default wth the KMV-Merton Model) provided the following expression for the implied probability of default, or EDF:

meanKMVMerton[t_]=Mean[kmvMERTONFirmValue[t]/.{(Subscript[v, 0])->assetValue, μ->defaultBoundary,σ->stdevReturns}];

KMVMertonEDF[μ_,σ_,debt_,t_]:=CDF[NormalDistribution[],(-(((Log[meanKMVMerton[t]])/debt)+(μ-0.5*(σ^2))*t)/(σ*Sqrt[t]))]

This model can now be provided directly to Manipulate to explore the credit rating assessment of the firm.

Manipulate[Plot[KMVMertonEDF[μ,σ,debt,t],{t,0,5},AxesLabel->{"Time Periods","Expected Default Frequency"}],{{μ,defaultBoundary,"Default Boundary"},10,120},{{σ,stdevReturns,"Standard Deviation of Returns"},2,50},{{debt,0.05,"Face Value of Firm's Debt"},0.01,0.1},SaveDefinitions->True]

One-Factor, Short-Interest-Rate Models

In all of the cases below, we give the SDE definition of the interest rate process, then its corresponding inbuilt Mathematica definition, show that they are the same, and plot multiple paths.

Merton's model (1973)

Merton’s model dr(t)==μdt+σ dw(t) is the Wiener process with nonzero drift μ and diffusion σ, and nonzero initial value (Subscript[r, 0]):

ItoProcess[WienerProcess[μ,σ]]

Use differential notation dr(t)==μ dt+σ dw(t) with initial condition r(0)=(Subscript[r, 0]) to define the same process:

ItoProcess[\[DifferentialD]r[t]==μ \[DifferentialD]t+σ \[DifferentialD]w[t],r[t],{r,(Subscript[r, 0])},t,w\[Distributed]WienerProcess[]]

This is equivalent to a simple Ito process with initial conditions:

ItoProcess[{μ,σ},{r,(Subscript[r, 0])},t]

merton[μ_,σ_,r0_]=ItoProcess[{μ,σ,r[t]},{r,r0},t]

ListLinePlot[#,AxesLabel->{"T","R(t)"},PlotLabel->"Merton Interest Rate Process"]& @RandomFunction[merton[0.1,0.05,0.02],{0.,1.,0.01},10,Method->"Milstein"]

Vasicek model (1977)

Vasicek’s model dr(t)==θ (μ-r(t))dt+σ dw(t) is the Ornstein–Uhlenbeck process with nonzero long-term mean μ and volatility σ, mean reversion rate θ, and nonzero initial value (Subscript[r, 0]):

ItoProcess[\[DifferentialD]r[t]==θ(μ -r[t])\[DifferentialD]t+σ \[DifferentialD]w[t],r[t],{r,(Subscript[r, 0])},t,w\[Distributed]WienerProcess[]]

ItoProcess[OrnsteinUhlenbeckProcess[μ,σ,θ,(Subscript[r, 0])]]

vasicek[μ_,σ_,θ_,r0_]=OrnsteinUhlenbeckProcess[μ,σ,θ,r0]

ListLinePlot[#,AxesLabel->{"T","R(t)"},PlotLabel->"Vasicek Interest Rate Process"]& @RandomFunction[vasicek[0.1,0.01,1.5,0.12],{0.,1.,0.01},10]

The Ornstein–Uhlenbeck process has a well-defined StationaryDistribution:

StationaryDistribution[OrnsteinUhlenbeckProcess[μ,σ,θ]]

It is the same for the process with any starting value:

StationaryDistribution[OrnsteinUhlenbeckProcess[μ,σ,θ,a]]

Rendleman–Bartter model (1980)

The Rendleman–Bartter model (1980) dr(t)==θ r(t)dt+σ r(t) dw(t) is the geometric Brownian motion process with nonzero drift θ, volatility σ, and nonzero initial value (Subscript[r, 0]):

ItoProcess[\[DifferentialD]r[t]==θ r[t] \[DifferentialD]t+σ r[t] \[DifferentialD]w[t],r[t],{r,(Subscript[r, 0])},t,w\[Distributed]WienerProcess[]]

ItoProcess[GeometricBrownianMotionProcess[θ,σ,(Subscript[r, 0])]]

rendlemanBarter[θ_,σ_,r0_]=GeometricBrownianMotionProcess[θ,σ,r0]

ListLinePlot[#,AxesLabel->{"T","R(t)"},PlotLabel->"Rendleman-Bartter Interest Rate Process"]& @RandomFunction[rendlemanBarter[0.2,0.5,0.1],{0.,1.,0.01},10]

Cox–Ingersoll–Ross model (1985)

The Cox–Ingersoll–Ross model dr(t)==θ (μ-x(t))dt+σ (Sqrt[x(t)]) dw(t) with long-term mean μ, volatility σ, speed of adjustment θ, and initial condition (Subscript[r, 0]):

ItoProcess[\[DifferentialD]r[t]==θ(μ -r[t]) \[DifferentialD]t+σ (Sqrt[r[t]]) \[DifferentialD]w[t],r[t],{r,(Subscript[r, 0])},t,w\[Distributed]WienerProcess[]]

CoxIngersollRossProcess[μ,σ,θ,(Subscript[r, 0])] represents a Cox–Ingersoll–Ross process with long-term mean μ, volatility σ, speed of adjustment θ, and initial condition (Subscript[r, 0]):

ItoProcess[CoxIngersollRossProcess[μ,σ,θ,(Subscript[r, 0])]]

ListLinePlot[#,AxesLabel->{"T","R(t)"},PlotLabel->"Cox–Ingersoll–Ross Interest Rate Process"]& @ParallelTable[RandomFunction[CoxIngersollRossProcess[0.03,0.4,0.3,0.04],{0.,1.,0.01}],{3}]