# Share experiences, preferably the surprising ones, with using ItoProcess

Can people please share their experiences, preferably the surprising ones, with using ItoProcess?

I am a big fan of ItoProcess and have already used it for several finance-related tasks, though I suspect that I'm just seeing the tip of the iceberg. I'll kick things off by supplying one answer about a bug that must be worked around before you can use ItoProcess with a correlation matrix. Wolfram Support have responded to my bug report by saying it will be fixed 'in a future release of Mathematica'.

• So could you please share one of your "several finance-related" experiences with us? – Rod Jun 15 '13 at 21:15
• Can you answer this question too? – cormullion Jun 15 '13 at 21:20
• Interesting, but questions posted here should generally be more concrete and less open ended than what you posted here. This is not a discussion forum. Could you try to rewrite your question and give it somewhat more focus? – Sjoerd C. de Vries Jun 15 '13 at 21:53
• Since this is going to seem like a motley collection of reports on ItoProcess[], I'm setting this CW for the time being. – J. M. will be back soon Jun 16 '13 at 1:48

## 3 Answers

Although I have found ItoProcess an excellent tool in nearly all cases, in those situations where a correlation matrix is supplied to it ItoProcess incorrectly uses the Cholesky decomposition of the correlation matrix when it should actually be using the transpose of the Cholesky. Note that I discovered this bug back in December 2012 when I was looking at the "Heston Model" example at this URL: http://reference.wolfram.com/mathematica/example/HestonModel.html

Before addressing the ItoProcess bug, let's revisit the problem of transforming two normally distributed random variables, y1 & y2, with correlation cc into two equivalent uncorrelated random variables, x1 & x2 via the use of the Cholesky decomposition…

Clear[ρ, x1, x2];
cc := {{1, ρ}, {ρ, 1}};
ccChol = Refine[CholeskyDecomposition[cc], -1 <= ρ <= 1]

Out := {{1, ρ}, {0, Sqrt[1 - ρ^2]}}


To turn this into a concrete example, prepare a sample of data representing y1 & y2…

n01 = {y1,y2} = RandomVariate[NormalDistribution[0, 1], {10^5, 2}];


We can see these are practically uncorrelated because Correlation[n01] yields something close to the Identity matrix.

OK, now we use the transpose of the Cholesky matrix to generate a sample of data, corset, with correlation, say, 0.75, which can be confirmed by using Correlation[corrset]…

ccChol75 = ccChol /. ρ -> .75;
corrset = Map[Transpose[ccChol75].# &, n01];


The symbolic version of the above is…

{y1, y2} = Transpose[ccChol].{x1, x2};


which results in the two transformed equations for y1 & y2 which reflect the required correlation…

{y1, y2} = {x1, x1 ρ + x2 Sqrt[1 - ρ^2]}


Now, we encapsulate the above with ItoProcess (using Refine to constrain the correlation to a sensible range) and demonstrate the bug with the use of the Correlation function…

ip = Refine[
ItoProcess[{{0, 0}, IdentityMatrix[2]}, {{w1, w2}, {0, 0}}, t, {{1, ρ}, {ρ, 1}}],
ρ - 1 <= ρ <= 1];
Correlation[ip[1],1,2]


ρ / Sqrt[1 + ρ^2]

This is clearly wrong - it should be simply ρ.

To explore the bug further, take a look at the structure of "ip"…

ItoProcess[{{0, 0}, {{1, ρ}, {0, Sqrt[1 - ρ^2]}}, {w1[t], w2[t]}}, {{w1, w2}, {0, 0}}, {t, 0}]

The clue to the problem is that the diffusion component of "ip" is the Cholesky matrix, whereas it should be its transpose. You can fix it by manually replacing the diffusion component in "ip" with the following piece of workaround code...

ItoProcessDougs[a_, b_, c_, Σ_?MatrixQ] := Module[
{ccChol, ip, rhoBounds, dum, nr, nc},
{nr, nc} = Dimensions[Σ];
rhoBounds =
dum[Flatten[Table[-1 <= Σ[[r, c]] <= 1, {r, nr}, {c, nc}]]] /. List -> Sequence;
rhoBounds = rhoBounds /. dum -> And;
ip = Refine[ItoProcess[a, b, c, Σ], rhoBounds];
ip[[1, 2]] = Transpose[ip[[1, 2]]];(*repair the ItoProcess*)
Return[ip]
];


Here is some more code which takes the ItoProcess and reconstitutes it back into the classic SDE format, which I find easier to work with if doing something like deriving, say, a stochastic-local volatility variant of the Heston model from first principles for use in, say, NDSolve (which can do quite a good job of numerically solving the resulting option pricing PDE, btw)...

ItoDriftCoefficient[pr_ItoProcess] := pr[[1,1]];
ItoDiffusionCOefficient[pr_ItoProcess] := pr[1,2]];

ToSDEs[proc_ItoProcess, {t_,ws_}] := Module[
{drifts = ItoDriftCoefficient[proc],
diffs = ItoDiffusionCoefficient[proc],
nEqn, nRV
},
nEqn = Length@drifts;
nRV = Length@First@diffs;
Table[drifts[[i]] dt + Sum[diffs[[i, j]] dws[[j]][t], {j,nRV}], {i,nEqn}]
];

ToSDEs[proc_ItoProcess] := Module[
{ws,nRV=Length@First@ItoDiffusionCoefficient[proc]},
ws=Table[w[j],{j, nRV}];
ToSDEs[proc, {t, ws}]
];


Applying the above to (buggy) "ip" (with the dt and dws terms appropriately differential) results in the incorrect output…

ToSDEs[ip, {t, {w1, w2}}]

Out := {dw1 + ρ dw2, Sqrt[1 - ρ^2] dw2}


Applying it to the "fixed" ItoProcess results in the correct output…

ipFixed =
ItoProcessDougs[{{0, 0}, {{1, 0}, {0, 1}}}, {{w1, w2}, {0, 0}}, t, {{1, ρ}, {ρ, 1}}]
Correlation[ipFixed[1], 1, 2]
ToSDEs[ipFixed, {t, {w1, w2}}]


ItoProcess[{{0, 0}, {{1, 0}, {ρ, Sqrt[1 - ρ^2]}}, {w1[t], w2[t]}}, {{w1, w2}, {0, 0}}, {t, 0}]
ρ
{dw1[t], ρ dw1[t] + Sqrt[1 - ρ^2] dw2[t]}

• This is a very interesting discovery! I would like to know if Wolfram has confirmed this bug for the 2-case ItoProcess[]. However, a localized bug does not invalidate the whole ItoProcess[] function in Mathematica, once for almost all the other cases the function works as intended. Furthermore, to find a localized bug is not a "surprising experience"... I was expecting something to help people understand the true utility of ItoProcess[] and practical applications for the function as well (especially in the Finance field !!!)... Anyway, great job! – Rod Jun 15 '13 at 22:27
• Thanks for the feedback, Rod. I plan to post a couple of things where ItoProcess actually does a fantastic job, but first things first. BTW, Wolfram have confirmed the bug. – doug fromoz Jun 15 '13 at 22:35
• And it seems this bug has been fixed in Mathematica 11. – Trad Dog Sep 16 '16 at 8:43

Derivation of the Heston model from first principles using ItoProcess[ ]

It's pleasantly surprising that something as subtle and useful as the Heston PDE can be easily derived via ItoProcess, once a few hurdles are overcome along the way. It's useful because, although the std Heston derivation is retold intuitively and in detail by Wilmott, by Shaw, by Gatheral, and many others, once you have coded the basic derivation machinery in Mathematica you can use it safely and robustly to do derivations of your own proprietary PDE models with, say, "stochastic-local" volatility. And even before you start the PDE derivation you can do some very powerful explorations of the system of SDEs by, say, using the ItoProcess directly in the calculation option prices via monte-carlo, or by exploring symbolic correlation formulae of the components. And then once the derivation is complete, the resulting PDE can be used directly in NDSolve[ ] to calculate prices via, say, finite-difference or method-of-lines methods, or the PDE can be (after reformatting) pasted directly into a tool like SciFinance's SciPDE.

Initial reference: http://www.wolfram.com/mathematica/new-in-9/time-series-and-stochastic-differential-equations/heston-model.html (note that, as-at this post, it uses the v9.0 version of ItoProcess, and therefore the graphics at that URL are subtly incorrect due to the correlation matrix bug).

But first let me take this opportunity to acknowledge Mark Fisher. Before ItoProcess[ ], for many years I used his Mathematica package which can be found these days at http://www.markfisher.net/~mefisher/mma/ItosLemma.m. In some ways Mark's package is still preferable to ItoProcess[ ] because it allows you to work with correlated wieners whereas ItoProcess[ ], when it can, reduces the system of equations down to orthogonal processes. That causes some initial discomfort when you first try and match some canonical examples, say Wilmott or Shaw or Gatheral. Enough waffle…

Following the authors noted above, we start the derivation by defining the ItoProcess for two correlated Wieners. Note that this applies the function defined in my earlier post, ItoProcessDougs[ ], which fixes a bug in the built-in ItoProcess for correlated processes...


cW[ρ_] := ItoProcessDougs[{{0, 0}, IdentityMatrix2}, {{w1, w2}, {0, 0}}, t,
{{1,ρ}, {ρ, 1}}];
cWproc = cW[ρ]
Correlation[cWproc1, 1, 2] ;
Out := ItoProcess[{{0, 0}, {{1, 0}, {ρ, Sqrt[1 -ρ^2]}},
{w1[t], w2[t]}}, {{w1, w2}, {0, 0}}, {t, 0}]
ρ


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


hm = ItoProcess[{ ds[t] == (μ - q) s[t]dt + Sqrt[V[t]] s[t]dws[t],
dV[t] ==κ (θ - V[t])dt +ξ Sqrt[V[t]]dwν[t]},
{s[t], V[t]}, {{s, V}, {s0, V0}}, t,
{ws, wν}[Distributed] cW[ρ]]
ToSDEs[hm, {t, {ws, W}}](*where ws and W are uncorrelated Wieners*)

Out := ItoProcess[{{-q s[t] +μ s[t],θκ -κ V[t]},
{{s[t] Sqrt[V[t]], 0}, {ξρ Sqrt[V[t]],ξ Sqrt[1 -ρ^2] Sqrt[V[t]]}},
{s[t], V[t]}}, {{s, V}, {s0, V0}}, {t, 0}]

{dt (-q s[t] +μ s[t]) +dws[t] s[t] Sqrt[V[t]],
ξ Sqrt[1 -ρ^2]dW[t] Sqrt[V[t]] +ξρdws[t] Sqrt[V[t]] +dt (θκ -κ V[t])}


Try out the SDEs in simulations and a monte-carlo option pricing

Simulate the model using the stochastic Runge-Kutta scheme :


hestClassic = {μ -> 0.1, q -> 0,κ -> 2,θ -> 1,ξ -> 1/2,ρ -> -1/3, s0 -> 25,
V0 -> 1.25};
τ = .25;
hestParams = {r ->μ, vκ ->κ, vL -> Sqrt[θ], vv ->ξ, vρ ->ρ, s -> s0,
v0 -> Sqrt[V0]} /. hestClassic;
td = BlockRandom[SeedRandom[1988];
RandomFunction[hm /. hestClassic, {0,τ, 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]]



optionParams = {p -> 1, k -> 3*s0, t ->τ} /. hestClassic;
optionParams = Union[optionParams, hestParams, hestClassic];

td = RandomFunction[hm /. hestClassic, {0,τ, 0.005}, 10^4,
Method -> "StochasticRungeKutta"];
pricePaths = td["PathComponent", 1];
Δ = pricePaths[τ];(*terminal distribution*)
optPayoffs = Map[Max[0, p (# - k)] /. optionParams &,
DistributionDomain[Δ]] /. optionParams;

(*calculate option price & std error*)
{simHes, simHesSE} = E^(-r t) {Mean[optPayoffs],
StandardDeviation[optPayoffs]/Sqrt[Length@optPayoffs]} /. optionParams

Out := {0.184978, 0.0235349}


Sorry, here I use one of my proprietary pricing functions, "Heston". But the interested reader can find useful algorithms at Shaw: http://www.mth.kcl.ac.uk/~shaww/web_page/papers/Heston09.pdf, and in Castagna's book, for example.


vvTiny = .00001;
priceHes = Heston[p, t, k, s, {v0, vv, vL, vκ, vρ}, q, r,λ] /. optionParams
priceHesBsm = Heston[p, t, k, s, {v0, vvTiny, v0, vκ, vρ}, q, r,λ] /. optionParams
priceBsm = BlackScholesMerton[p, t, k, s, v0, q, r] /. optionParams
Print[TableForm[{{"BSM", "SemiAn Hest", "Sim Hest", "Sim Hest SE"},
{priceBsm, priceHes, simHes, simHesSE}}]];

Out:=  0.172012
0.243079
0.242999


Derive the Heston PDE from the SDEs using BSM risk-neutral argument  (HELP!: I ran into trouble trying to render Mathematica's pd superscript format - can someone please provide a hint how to do it? I've left the mangled pd superscript code in this post just in case it's useful to someone...)

First, calculate the stochastic derivative of U(s,v,t) by taking the ItoProcess of the ItoProcess…  (Note: because ItoProcess reduces the processes to orthogonal rather then correlated wieners, we depart temporarily from the canonical derivations mentioned above)


Uproc = ItoProcess@ItoProcess[{ ds[t] == (μ - q) s[t]dt + Sqrt[v[t]] s[t]dws[t],
dv[t] ==κ (θ - v[t])dt +ξ Sqrt[v[t]]dwν[t]},
U[s[t], v[t], t], {{s, v}, {s0, v0}},      t,
{ws, wν}[Distributed] cW[ρ]];
dU = Last@ToSDEs[Uproc, {t, {ws, W}}];(*\W is a Wiener uncorrelated to ws*)
dU = Collect[dU,SuperscriptBox[U, TagBox[RowBox[{"(",
RowBox[{"", ",", "", ",", "_"}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t]]

(The following is a GIF image of the above mangled expression)

Now, (unnecessarily) use the fact that dwν[t] =ρdws[t] + Sqrt[1 -ρ^2]dW[t] to reconstitute dwν[t] - this is just to show the SDE using the traditionally employed correlated "volatility" Wiener.


dU /.ξ ρ dws[t] Sqrt[v[t]] +ξ Sqrt[1 -ρ^2]dW[t] Sqrt[v[t]] ->ξ Sqrt[v[t]]dwν[t]


Now calculate dS and dU1...


Sproc = ItoProcess@ItoProcess[{ ds[t] == (μ - q) s[t]dt + Sqrt[v[t]] s[t]dws[t],
dv[t] ==κ (θ - v[t])dt +ξ Sqrt[v[t]]dwν[t]},
s[t], {{s, v}, {s0, v0}}, t,
{ws, wν}[Distributed] cW[ρ]]

Out := ItoProcess[{{-q s[t] +μ s[t],θκ -κ v[t]},
{{s[t] Sqrt[v[t]], 0}, {ξρ Sqrt[v[t]],ξ Sqrt[1 -ρ^2] Sqrt[v[t]]}},
s[t]}, {{s, v}, {s0, v0}}, {t, 0}]


dS = First@ToSDEs[Sproc, {t, {ws, W}}](*\W is a Wiener uncorrelated to ws*)

Out:= dt (-q s[t] +μ s[t]) +dws[t] s[t] Sqrt[v[t]]


U1proc = ItoProcess@ItoProcess[{ ds[t] == (μ - q) s[t]dt + Sqrt[v[t]] s[t]dws[t],
dv[t] ==κ (θ - v[t])dt +ξ Sqrt[v[t]]dwν[t]},
U1[s[t], v[t], t], {{s, v}, {s0, v0}}, t,
{ws, wν}[Distributed] cW[ρ]];
dU1 = Last@ToSDEs[U1proc, {t, {ws, W}}];(W is a Wiener uncorrelated to ws)
dU1 = Collect[dU1,SuperscriptBox[U1, TagBox[RowBox[{"(",
RowBox[{"", ",", "", ",", "_"}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t]]


Ito derivative, dΠ, of the derivative value, U[s[t],t], hedged with Δ of stock, s[t] and Δ1 of a 2nd instrument, U1[s[t],t] , form a portfolio of one "short" derivative plus a "long" position of Δ units of stock, plus a "long" position of Δ2 units of another derivative with exposure to the stock and its volatility


Clear[Δ, V, s, v, r,σ,μ,κ,ξ, t, s0, w1, w2];
Π = U[s[t], v[t], t] -Δ s[t] -Δ1 U1[s[t], v[t], t];
dΠ = dU -Δ (dS + s[t] qdt) -Δ1 dU1 // FullSimplify


Extract the two risky components of dΠ, solve for Δ & Δ1, then substitute back into dΠ


eqnSrisk = dΠ /. {dt -> 0,dW[t] -> 0}
eqnWrisk = dΠ /. {dt -> 0,dws[t] -> 0}


find the values of Δ & Δ1 which completely hedges the stochasitic risk of the derivative


solnHedge = First@FullSimplify@Solve[eqnSrisk == 0 && eqnWrisk == 0, {Δ,Δ1}]


substitute Δ & Δ1 back into dΠ to make it "riskless" (by removing all the dw terms)


dΠ = FullSimplify[dΠ /. solnHedge]


Apply the risk-neutrality argument.  Now, because we've shown that the hedged portfolio is riskless, apply the Black-Scholes risk-neutrality argument and derive the Black-Scholes PDE by recognising that the return on the hedged portfolio must be the risk-free rate.


mma = r (Π /. solnHedge); (* risk-free money market account*)
lhs = FullSimplify[dΠ/dt - mma]; (*i.e. left hand side of dΠ-mma*)
pde = lhs == 0


Collect the U and U1 terms to observe that each "side" is equivalent to some function of S, v & t


Clear[v];
pde = Collect[pde, {SuperscriptBox[U, TagBox[RowBox[{"(",
RowBox[{"", ",", "", ",", ""}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t],
SuperscriptBox[U1, TagBox[RowBox[{"(",
RowBox[{"", ",", "", ",", ""}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t]}]


From the above, observe that collecting terms in U to one side of the equation, and U1 terms to the other, results in two equal expressions that have identical coefficients. Hence, we deduce that each expression is equal to the same equivalent function of S, v & t. So we arbitrarily define that equivalent function in terms of the expression containing only terms in the "extra" asset price, U1, and its greeks.

So, isolate the ∂U/∂v term from the above as the "equivalent function"...


Clear[equivFn];
equivFn = First[pde /.SuperscriptBox[U, TagBox[RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t] -> 1 /.SuperscriptBox[U,
TagBox[RowBox[{"(",RowBox[{"", ",", "", ",", "_"}],
")"}],Derivative],MultilineFunction->None][s[t], v[t], t] -> 0
/.SuperscriptBox[U, TagBox[RowBox[{"(", RowBox[{"0", ",", "0", ",", "0"}],
")"}],Derivative],MultilineFunction->None][s[t], v[t], t] -> 0]


Finally, replace the expression representing the "equivalent function" with an allowable functional form that includes what's called the "market price of volatility risk", λ ,to produce the final Heston PDE. λ is called the market price of volatility risk because it tells us how much of the expected return of U is explained by the risk (i.e. standard deviation) of v in the Capital Asset Pricing Model framework.  The "equivalent function" is chosen to also include the coefficient of the first derivative of v in the dΠ equation, before hedging. This completes the derivation of the Heston PDE.


Clear[α,β,λ];
pde = Collect[pde, {SuperscriptBox[U, TagBox[RowBox[{"(",
RowBox[{"", ",", "", ",", ""}], ")"}],Derivative],
MultilineFunction->None][s[t], v[t], t],
SuperscriptBox[U1, TagBox[RowBox[{"(",RowBox[{"", ",", "", ",", ""}],
")"}],Derivative],
MultilineFunction->None][s[t], v[t], t]}]
/. {equivFn -> (κ (θ - v[t]) -λ[s[t], v[t], t])} // FullSimplify


Reformat the Heston PDE to more conventional form using a (slightly tweaked version of a) tool to be found on the web...


Clear[S];
FormatPDiff["Classic"];
pde /. {s[t] -> S, v[t] -> v}
FormatPDiff[];

Out := 2 (r U[S, v, t] + (v κ -θ κ +λ[S, v, t])∂U/∂v + (q - r) S ∂U/∂S)
==   2 ∂U/∂t + v (ξ^2 ∂^2U/∂v^2 + S (2ξ ρ ∂^2U/∂S∂v + S ∂^2U/∂S^2))


In a followup post I intend to explore a couple of things touched on above - maybe a comparison of the ItoProcess "Stochastic Runge-Kutta" vs. Andersen's "Quadratic-Exponential", and maybe using the derived PDE in NDSolve[ ]. Anyone else wants to jump in, too, please be my guest.

• I just took a look at this in an IE browser - for some reason the <pre><code> this=some_code; </pre></code> formatting is not rendering properly, whereas on my Mac (and iPhone) it does - any hints? – doug fromoz Jun 26 '13 at 23:03
• Probably not what you were after but the pd display issue can be worked around by using the full form, e.g. dU = Collect[dU, Derivative[_, _, _][U][s[t], v[t], t]]. By the way, since 2013, it seems things have improve and the Collects are not needed to obtain your results. – Trad Dog Sep 16 '16 at 10:38
• Very interesting example of the use of ItoProcess. Thanks for the extensive derivation. – Trad Dog Sep 16 '16 at 10:39
• By the way, with the version of FormatPDiff, I do not get as nearly a result as yours. Would you mind sharing the "tweak". – Trad Dog Sep 16 '16 at 10:41

We can take the Heston equation, 'pde', from the previous post and do some neat calculations with it, via NDSolve[ ]. Now, I don't normally use NDSolve[ ] for this kind of stuff, so if there's a better way of doing things, please don't hesitate to let us know... And sorry, I still haven't gotten the hang of entering raw Mathematica code into this wiki environment, so I'm afraid I am going to use GIFs again.

Firstly, for convenience, assign some of Heston's 'classic' greek-style variables to ascii-based names (so that, eg, VBA & C++ code is somewhat interchangeable). And again I use some proprietary software to calculate the Heston and Black-Scholes-Merton option prices for the given parameters, but please see the references mentioned above in the previous post to download equivalent models. The last two prices of about 24.3 cents are the Heston price with negligible stochastic vol and std BSM respectively (i.e. essential equivalent), whereas the 3rd-last value of about 17.2 cents is the heston call price in the presence of significant stochastic vol...

So essentially our goal is to use NDSolve[ ] to replicate that heston price of about 17.2 cents.

Because we are going to be using numerical methods, we need to have an idea of the space in which we are working. That is, the min and max asset prices and volatilities that we will be working within when calculating our solution. The following shows how to use ItoProcess[ ] to get an estimate of the "6 standard deviation" upper and lower bounds: SMin, SMax, VMin and VMax (the last two being variance rather than volatility). When taken with the time constraints of "now" (ie. t=0) and the option expiry, these bounds effectively put the Ito process in a 3D box (i.e. denote the boundary conditions) which is then solvable with NDSolve[ ]...

Now that we have the pde (which we've earlier mapped to "pde2" so we can use convenient symbols) constrained within reasonable bounds, it is then fairly trivial to put it into NDSolve[ ] and get a Heston option price. I'm pretty much a novice at using NDSolve[ ], so if you know of better, more accurate or efficient ways of doing this, please don't hesitate to jump in and let us know...

OK, there are some warnings - PLEASE HELP ME TO MAKE THEM GO AWAY! :) And the graphical surface of option value looks ghastly - help with that would be welcomed, too. But notwithstanding the awkwardness of the solution method, it gives a value close to the expected value of 17.2 cents (which is significantly different to the BSM value of 24.3 cents), so something must be working ok. You can see that I've also tried "method of lines" of the pseudospectral flavour, and I've tried a few other NDSolve options, though I've found (so far) that the extra computation time is not worth the additional precision (IMHO).

Well, this nice bottle of red is bottoming, so I'll leave it at that. Cheers.