Message Boards Message Boards

20 Replies
9 Total Likes
View groups...
Share this post:

Solve the probability distribution for the Forward Kolmogorov.

Posted 10 years ago

Hello everybody,

I would be forever grateful if you can get me started on how to solve the following:

enter image description here

I understand how to go from the stochastic differential equation to the Forward Kolmogorov; but, I do not know how to solve for p(t,x). Furthermore, the textbooks that I read specify that at t=0, the probability distribution is the DIrac Delta function. I am assuming this implies that we integrate the system of equations somehow?!

Also important for me is to learn how to define the joint density (the solution) so I can use the Maximum Likelihood to fit a given set of data.

Thank you in advance! Again - any help would be greatly appreciated.


POSTED BY: Edvin Beqari
20 Replies
Posted 10 years ago

Mr. Lichtblau, thank you for your guidance. I am tracking 100% with your reply above. The poster above (Henrik) shed some light on the FT conversion of the unknown function - even though there is still an issue as I specified in my reply to him.

Nevertheless, I was able to convert my hand derivation with Mathematica and it scaled properly so I am hopeful we are almost there. I am sorry again for not being clear since the beginning.

POSTED BY: Edvin Beqari
Posted 10 years ago

Ok, thank you sir. I will try to show again what I am attempting to do.

Step 1: Define an Ito Process (Regular Brownian Motion in this case)

bm[\[Mu]_, \[Sigma]_, x0_] := 
     t] == \[Mu]* \[DifferentialD]t + \[Sigma] \[DifferentialD]w[t], 
  x[t], {x, x0}, t, w \[Distributed] WienerProcess[]]

Step 2: Extract the Transition Probability Density

PDF[bm[\[Mu], \[Sigma], x0][t], x]


E^(-((x - x0 - t \[Mu])^2/(2 t \[Sigma]^2)))/(Sqrt[2 \[Pi]] Sqrt[
 t \[Sigma]^2])

Step 3. Define a slightly more complicated Ito Process (Birth and Death diffusion in this case)

bd[\[Beta]_, \[Delta]_, x0_] := 
     t] == ((\[Beta] - \[Delta])*x)* \[DifferentialD]t + 
    Sqrt[(\[Beta] + \[Delta])*x]*\[DifferentialD]w[t], x[t], {x, x0}, 
  t, w \[Distributed] WienerProcess[]]

Step 4. Extract the Transition Probability Density

PDF[bd[\[Beta], \[Delta], x0][t], x]


"Parameter 0 at position 1 in \
NoncentralChiSquareDistribution[0,-((4\E^(t (\[Beta]-\[Delta]))\x0\(\
\[Beta]-\[Delta]))/((1-E^(t \
(\[Beta]+Times[<<2>>])))\(\[Beta]+\[Delta])))] is expected to be \

So, the above example illustrates that there are some interesting diffusion that the inversions Mathematica uses do not work. Therefore, there is a need for another method to extract the density function.

Another method involves - Deriving the Forward Kolmogorov for the specified diffusion. Then, reduce from a second order to a first order with Fourier Transforms - apply, the initial condition - and inverse back to get the function. - I would like to be able to do this in Mathematica. I do not have a code for this yet.*

So, this is the whole story in a nutshell.

Udo Krause above is on it - It just needs a little bit of modification perhaps to get it in the right form.

k3[\[Tau]_, \[Xi]_] := 
 1/Sqrt[2 \[Pi] \[Tau]] Exp[-((\[Xi] - \[Xi]0 + a f[\[Tau]])/
        b)^2/(2 \[Tau])]

D[k3[t, x], t] + a D[k3[t, x], x] - 
  b^2/2 D[k3[t, x], {x, 2}] // Simplify


-((a E^(-((x - \[Xi]0 + a f[t])^2/(
   2 b^2 t))) (x - \[Xi]0 + a f[t]) (1 + Derivative[1][f][t]))/(
 b^2 Sqrt[2 \[Pi]] t^(3/2)))

Then he solves for f' but - I don't follow anymore from that point.

Nonetheless, I would like to get the Fourier Transforms method right with Mathematica as it a very nice tool. I hope this clarifies the issue a little, and I apologize for not being very clear before.

POSTED BY: Edvin Beqari

This is a big improvement in the sense that it has actual code and a more clearly stated question.

Even with FT methods, I'm not sure you will get far unless you can figure out how to evaluate the transform in question. It seems to be difficult at least if tackled as an integral.

The error message from PDF[bd[\[Beta], \[Delta], x0][t], x] is due to there being no constant term inside the square root in bd[...]. If you add a smallish value (I use 1/1000 below) then that goes away.

bd[\[Beta]_, \[Delta]_, x0_] := 
     t] == ((\[Beta] - \[Delta])*x[t])*\[DifferentialD]t + 
    Sqrt[1/1000 + (\[Beta] + \[Delta])*x[t]]*\[DifferentialD]w[t], 
  x[t], {x, x0}, t, w \[Distributed] WienerProcess[]]

Alas, this still does not evaluate: PDF[bd[\[Beta], \[Delta], x0][t], x]. So I'm not sure it helps to make this change.

POSTED BY: Daniel Lichtblau

I would suggest posting only Mathematica code. Even if Mathematica cannot evaluate it, at least readers will know what you have tried and also have a starting point for perhaps going further with it. From the posts above I have some idea of what you have tried but it is far from a complete picture and it is not obvious what is wanted.

As it now stands, I see a bunch of images, one of which is a PDE that I will surmise you wish to solve. I also see images of classical Fourier transform formulas. It is quite unclear to me what it is you want, how the transforms fit into the PDE, and so forth. If your interest is in using FT methods to solve a PDE you should state that.

Also please note that an image of actual code is of significantly less value than the actual code itself in cut-and-paste form.

POSTED BY: Daniel Lichtblau
Posted 10 years ago

My ultimate goal, as I have expressed earlier in the thread, is to be able to solve the equation using Mathematica.

  • My first example (pictures of notes) is an example of something that Mathematica can solve.Hence I showed the math and the code how to obtain the same result with Mathematica.

  • My second example (picture of notes 2 ) is an example that Mathematica cannot solve using the same technique, so its not simply a rehash (Weiner Process vs Birth&Death Diffusion)

If I do not show the mathematics behind it - how can I write a meaningful code?

Nonetheless, If I am posting in the wrong group - I apologize. Please direct me in the right one, or move/ delete the thread.

POSTED BY: Edvin Beqari
Posted 10 years ago

I also want to know how to define the unknown function p(x,t) in Mathematica so I can take the FT of its derivative w.r.t to x and t? - Then solve it as a regular ODE with Dsolve including the initial condition.

  • All the examples in Wolfram Documentation are FTs of know functions.
POSTED BY: Edvin Beqari

Dear Edvin,

posting pictures of basic Fourier identities is not conducive to good questions. This forum is first of all about Wolfram Technologies. Trying to improve your Wolfram Language code and asking questions about your Wolfram Language code is the proper approach.

Also your new notes should really contain more than just a rehash of earlier ones.

Sincerely, Moderation Team

POSTED BY: Moderation Team
Posted 10 years ago

enter image description here

enter image description here

POSTED BY: Edvin Beqari

Dear Edvin,

I did not follow this thread very closely, but I would like to make a remark concerning Fourier transform:

There is no specific "Mathematica definition" of Fourier transform, but Mathematica follows - as standard setting - just the standard notation of classical physics. If one does not like this (e.g. myself) it can be changed - as it is well described in the documentation ("Details and Options").

I am referring now to the emphasized formulas above "Definitions & Properties Table": According to the definition given in the first line you need to change to the notation of "signal processing" by setting the option

FourierParameters -> {0, -2 \[Pi]}

Then amazing things can be done, e.g.:

In[1]:= FourierTransform[f''[x], x, k, FourierParameters -> {0, -2 Pi}]

Out[1]= -4 k^2 \[Pi]^2 FourierTransform[f[x], x, k]

In[2]:= InverseFourierTransform[-I  g'[k], k, x, 
 FourierParameters -> {0, -2 \[Pi]}]

Out[2]= -2 \[Pi] x InverseFourierTransform[g[k], k, x]

In[3]:= FourierTransform[DiracDelta[x - x0], x, k, 
 FourierParameters -> {0, -2 Pi}]

Out[3]= E^(-2 I k \[Pi] x0)

From this one can see that the third formula is not consistent with the definition (first line). Furthermore the second formula does not make sense anyway ...

Cheers Henrik

POSTED BY: Henrik Schachner
Posted 10 years ago

Dear Henrik, thank you so very much. The conversion was so helpful and I was able to convert my derivation # (7) to the p(x,t). Consequently, I fixed t, and integrated over [-Inf, inf ], and it gave me 1 - which makes me believe it is the right answer.

I am not completely in the free yet. As one may notice, in (2) and (3) - I am taking the FT of the derivative of x * f(x), not of f(x).

In[8]:= FourierTransform[x*f''[x], x, k, 
 FourierParameters -> {0, -2 Pi}]

Out[8]= -2 I (2 k \[Pi] FourierTransform[f[x], x, k] + 
   k^2 \[Pi] FourierTransform[I x f[x], x, k])

I do not understand the output, and furthermore the x * f(x) term did not go away. I was only able to do it by hand because I found property (c) online (more than one source). Also, Definition (a) and (b) are straight from the book, and (d) checks out with the Mathematica output.

Nevertheless, it was really good to learn the {0,- 2Pi} conversion - although, I still cannot perform the whole conversion with Mathematica.

Thanks again Henrik.

PS: I updated the Table above as I found an error with my algebra.

POSTED BY: Edvin Beqari
Posted 10 years ago

So, I am back at working on this and I suspect I am on the right track this time.

The question I raised earlier in the thread was why do we need the delta function as the initial condition? - One answer that made sense to me is the following: - When we solve an ordinary diff eq we end up with a constant for which we apply the initial condition and solve for.

  • In the Kolmogorov equation above, p(t,x) is a function of x and t, but we are taking a partial in respect to time. Therefore:
    • If we do the Fourier Transform of each term - and reduce the PDE to an ODE - since p(t,x) is a function of two variable we need a function instead of a constant as an initial condition (Recall the a = b(x) trick if I remember correctly from my diff eq class).
POSTED BY: Edvin Beqari
Posted 10 years ago

Mr. Krause, I was wondering if you were able to figure out anything from the book above. I will follow up the thread with the "Characteristic Method," but I have not had much time lately. Thank you for all your help though.

POSTED BY: Edvin Beqari
Posted 10 years ago

OK, so allow me to back up a little and give a little more background...

I am interested in diffusion - not necessarily the heat diffusion. However, from what I understand, all diffusion are modeled using the heat equation. Furthermore, an extension of the derivative driven diffusions are the stochastic differential equation diffusions where the diffusion coefficient is perturbed with a Wiener Process.

Here it is (same as the one in my first post):

simpleDiff[\[Alpha]_, \[Beta]_, x0_] := ItoProcess[\[DifferentialD]x[t] == \[Alpha] \[DifferentialD]t + \[Beta] \[DifferentialD]w[t], 
  x[t], {x, x0}, t, w \[Distributed] WienerProcess[]]

With some diffusions such as the one above - Mathematica can extract the joint density distribution - which is nice because it allows us to use the Maximum Likelihood which is the preferred method to fit the parameters.

PDF[simpleDiff[\[Alpha], \[Beta], x0][t], x]

enter image description here

Now let's say I want to model some bird population which depends on the birth and death rate of this population - but all we can observe is the population number whether it increases or decreases. So, the drift term would be (birth - death = mu), and it can be shown that such model can be defined as follows:

birdpop[\[Mu]_, \[Sigma]_, x0_] := ItoProcess[\[DifferentialD]x[t] == \[Mu]*x \[DifferentialD]t + 
    Sqrt[\[Sigma]*x] \[DifferentialD]w[t], x[t], {x, x0}, t, w \[Distributed] WienerProcess[]]

But if we try to get the joint density as in the first example - it just does not work. So, is the case for many other diffusion models that I have looked at.

PDF[birdpop[\[Mu], \[Sigma], x0][t], x]

This is where the need arises to derive the Forward Kolmogorov Equation of the SDE, then solve the partial differential equation with the initial condition as the Dirac Function. There is not much explanation in the books I have looked at why this is the case. I am just following their examples in attempt to understand. All I know is that I have to think of the DiracDelta as the limit of the "unknown distribution" as t->0.

Furthermore, since we have two independent variables in the PDE [ t and x ] - we have to use the Fourier Transform which converts the PDE in an ODE ... just like the steps in my second post. Actually, most complex models have to be solved using the characteristic methods, but for now I wanted to get through the easier one first.

So, this is basically the whole picture for now. And, by no means, am I implying that what is written above is wrong - I just do not know how to utilize it with what I am doing. I have received the "Generalized Functions in Mathematical Physics' book so I will try to read up on it as time permits since school keeps me pretty busy already.

POSTED BY: Edvin Beqari

Okay, only now I've seen your post Fitting Data to Stochastic Differential Equation (ITO Process) related to this one ... and - oh happy day - the referenced book by Dr. Allen is available on the net Allen Itô ... let's see how they describe things.

POSTED BY: Udo Krause

First, with k4 also a k5 with an arbitrary power of b in front is a solution:

  In[1]:= Clear[k5, b, a, m]
    k5[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] b^m \[Tau]] Exp[- ((\[Xi] - \[Xi]0 - a \[Tau])/b)^2/(2 \[Tau])]

  In[3]:= D[k5[t, x], t] + a D[k5[t, x], x] - b^2/2 D[k5[t, x], {x, 2}] // Simplify
  Out[3]= 0

  In[4]:= D[k5[t, x - y], t] + a D[k5[t, x - y], x] - b^2/2 D[k5[t, x - y], {x, 2}] // Simplify
  Out[4]= 0

now, what is the solution k[t,x] for a general initial condition IC (see V. S. Vladimirov for specifications)?

k[t,x] = Integrate[k5[t,x-y] IC[y], dy]
D[k[t, x], t] + a D[k[t, x], x] - b^2/2 D[k[t, x], {x, 2}] = Integrate[D[k5[t, x-y], t] + a D[k5[t, x-y], x] - b^2/2 D[k5[t, x-y], {x, 2}] IC[y], dy]
                                                           = Integrate[0 IC[y], dy] = 0 for t > 0

for t = 0 there is

 Integrate[D[k5[t, x-y], t] + a D[k5[t, x-y], x] - b^2/2 D[k5[t, x-y], {x, 2}] IC[y] dy]
 = Integrate[DiracDeltat[x - y] IC[y], dy]
 = IC[x]

initial condition fulfilled.

What's your point with the DiracDelta?

The important property is

  In[11]:= Integrate[f[x - y] DiracDelta[y], {y, -Infinity, Infinity}]
  Out[11]= f[x]

  In[12]:= Integrate[DiracDelta[x - y] f[y], {y, -Infinity, Infinity}]
  Out[12]= ConditionalExpression[f[x], x \[Element] Reals]

You can not use DiracDelta like an orderly function. That's why it belongs to the generalized functions. Dirac invented it by his needs. Von Neumann said "There is no need for a delta function." Then Schwartz won a Fields Medal for making it intelligible to other people than theoretical physicists. So it can be used safely since then. The mean feature is to have some function space where the generalized functions can act on. All identities have to be written with this function space. Strictly speaking it's non-sense to say

operator[f][x] == DiracDelta[x]

that's slang and an abbreviation for the correct thing to say, which is

Integrate[operator[f][x] g[x], x] == Integrate[DiracDelta[x] g[x], x]

where the g[x] are from the above mentioned function space. It's really functional analysis. In a way you should become familiar with it in order not to run in problems people had before Schwartz worked it out.

POSTED BY: Udo Krause
Posted 10 years ago

Dear Udo,

Thank you so much for your reply again! You are awesome.

So to my understanding, at the last part - you have shown that the equation in my first post is indeed the solution to the PDE with the delta initial condition.

However, the solution (joint density in (t,x)) is unknown and what we are solving for. To demonstrate my point let me purposely do it wrong below. I am going to slightly change the notation to match the steps in my second post.

First lets define deltafunction as in k0 :

f[t_, x_] := 1/Sqrt[2 \[Pi] t] Exp[-(x - x0)^2/(2 *t)]

Take the Fourier transform and suppress the output because it pretty long (perhaps an indication that is very wrong) :

fhat[t_, s_] := FourierTransform[f[t, x], t, s]

Then lets try to use DSolve to get the solution in transform space: Note: I am trying to use the same steps as above.

DSolve[{D[uhat[s], {x, 1}] == (a^2) D[uhat[s], {x, 2}], 
   uhat[0] == fhat[t, s]}, uhat, s] // Simplify

And, as one may guess, we do not get anything that makes sense - and most importantly if I input a zero in fhat[0,s] I get the infinite expression which emphasizes the problem on how to define the DiracDelta. Perhaps, your reply already has the answer but I am unable to see it. I am sorry.

POSTED BY: Edvin Beqari
Posted 10 years ago

I think I found something useful in some course notes online, although - this is a simpler example of a pure diffusion. I retyped the equations and I will try to mimic the logic in Mathematica sometime this coming weekend.

Perhaps, If I replace the initial condition f(x) with DiracDelta[x-x0] it may work.

One thing I have noticed with DiracDelta function in Mathematica is that it does not output in the distribution form. This could be a problem.

If someone is interested to give this a head start -> more than welcome. This stuff is fun.

enter image description here

POSTED BY: Edvin Beqari

I'm afraid there is no shorter learning path than the lecture of Vladimirov's book.

I'm afraid I'm working in greyed software producing chain gangs too long already ...

The initial condition u[0,x]=DiracDelta[x-x0] means that in the limit t->0 the solution becomes a delta function. The delta function can itself be defined as a limit. Why not take an appropriate limit definition and fit it into the solution?

Start with k0

Clear[k0, b]
k0[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] \[Tau]] Exp[-(\[Xi] - \[Xi]0)^2/(2 \[Tau])]

In[28]:= D[k0[t, x], t] - b^2/2 D[k0[t, x], {x, 2}] // Simplify
Out[28]= ((-1 + b^2) E^(-((x - \[Xi]0)^2/(2 t))) (t - (x - \[Xi]0)^2))/(2 Sqrt[2 \[Pi]] t^(5/2))

there is a problem with the b, check the coefficient needed

Clear[k1, b, n]
k1[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] \[Tau]] Exp[- (\[Xi] - \[Xi]0)^2 b^n/(2 \[Tau])]

In[19]:= D[k1[t, x], t] - b^2/2 D[k1[t, x], {x, 2}] // Simplify
Out[19]= -(((-1 + b^(2 + n)) E^(-((b^n (x - \[Xi]0)^2)/(2 t))) (-t + b^n (x - \[Xi]0)^2))/(2 Sqrt[2 \[Pi]] t^(5/2)))

this gives 2 + n == 0. Now take the convection term a != 0 into consideration

In[31]:= Clear[k2, b, a]
k2[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] \[Tau]] Exp[- ((\[Xi] - \[Xi]0)/b)^2/(2 \[Tau])]

In[33]:= D[k2[t, x], t] + a D[k2[t, x], x] - b^2/2 D[k2[t, x], {x, 2}] // Simplify
Out[33]= (a E^(-((x - \[Xi]0)^2/(2 b^2 t))) (-x + \[Xi]0))/(b^2 Sqrt[2 \[Pi]] t^(3/2)) 

that did not work, of course. How to modify the evolving guess k3? For a = 0 the function k2 should result. Any term a f[?] would do. Should f depend on both - time and space - variables or on time variable alone? There is already a dependency on space variable in the term - just do something simple:

Clear[k3, b, a, f]
k3[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] \[Tau]] Exp[- ((\[Xi] - \[Xi]0 + a f[\[Tau]])/ b)^2/(2 \[Tau])]

In[77]:= D[k3[t, x], t] + a D[k3[t, x], x] - b^2/2 D[k3[t, x], {x, 2}] // Simplify
Out[77]= -((a E^(-((x - \[Xi]0 + a f[t])^2/(2 b^2 t))) (x - \[Xi]0 + a f[t]) (1 + Derivative[1][f][t]))/(b^2 Sqrt[2 \[Pi]] t^(3/2)))

okay, we are done if one could solve 1 + f'[t] == 0 and that's possible :-)

In[81]:=DSolve[D[f[t], t] == -1 , f , t]
Out[81]= {{f -> Function[{t}, -t + C[1]]}}


In[82]:= Clear[k4, b, a]
k4[\[Tau]_, \[Xi]_] := 1/Sqrt[2 \[Pi] \[Tau]] Exp[- ((\[Xi] - \[Xi]0 - a \[Tau])/b)^2/(2 \[Tau])]

In[84]:= D[k4[t, x], t] + a D[k4[t, x], x] - b^2/2 D[k4[t, x], {x, 2}] // Simplify
Out[84]= 0

the fundamental solution k4 has flippered out of the initial condition.

P.S. In eqn. (10) convolution integrates over y, not x.

POSTED BY: Udo Krause
Posted 10 years ago

Thank you very much for the reply. I followed your advise and purchased the book from Amazon. I will update the thread if I figure out something more.

POSTED BY: Edvin Beqari

This Forward Kolmogorov Equation is in physics known as heat conduction equation. The DiracDelta as initial condition lets one find a so-called fundamental solution. It's called fundamental solution because one gets a solution for other initial conditions - roughly speaking - by convolution of the fundamental solution with the current initial condition (with other words, you must not solve over and over again for different initial conditions). All this is described e.g. in Generalized Functions in Mathematical Physics by V. S. Vladimirov.

I'm afraid there is no shorter learning path than the lecture of Vladimirov's book.

POSTED BY: Udo Krause
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract