x - Noppa

Lectures 11: Probability Density solves
Fokker-Plank Equation
Topics:
1. Fokker-Plank equation
2. Interpretation of time-varying probabilities
3. Wealth characterisation (portfolio problem)
4. Transition probability. ODE instead of Fokker-Plank equation
5. Computing of probabilities with Comsol
Literature:
 Liptser, Shiryaev. Statistics of random processes. Part II:
Applications. 2000. Chap. 9
 Lecture notes
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
1
Introduction
A nonlinear stochastic process
dz  a (t , z )dt  b(t , z )dW
can be characterised with the probability density p ( t , x ) that satisfies the
Fokker-Plank Equation, also called as Kolmogorov Forward Equation (KFE)
p (t , x) 1  2 2



b
(
t
,
x
)
p
(
t
,
x
)
 a(t , x) p(t , x) 

2 
t
x
2 x
p (0, x)  p0 ( x)
Our aim is to understand and apply this equation in process characterization.
It is a relatively simple equation for computing of the unconditional probability
density, mean and variance for any nonlinear process z (t ) .
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
2
Characterization Method (better Method than Ito’s Rule
This is a more general method than stochastic process characterization by Ito’s rule.
In investment problem we wonder: How bad or good our wealth can go?
This is even better, if we know the full range of possible wealth values in between
some minimal and maximal limits. Our specific interest is to know: How far apart
the trajectories of wealth can go at a given time?
1.5
Hence, we are interested in the
dispersion of trajectories
1
0.5
We solve this problem in two cases:
1. Wiener process (random walk) in this Fig
0
-0.5
-1
2. General process later
-1.5
-2
0
0.2
0.4
0.6
Time, s
0.8
Time
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
3
1
Wiener Process: Dispersion of Trajectories
Although the mean of a Wiener process is zero, its variance rises
progressively since the Wiener process trajectories disperse rapidly
dz  dW
W (0)  0
The probability density as a measure that tells: how rapidly. This density can be
found from the Fokker-Plank equation
p (t , x) 1  2 p (t , x)

t
2 x 2
 , x  0
 Initial density
p (0, x)   ( x)  
0, x  0
x represents a range xi  xi  xi 1 where a possible
trajectory can be (between gridlines) at time t:
xi 
The probability density p (t , x ) is a set function f ( xi ) );
it characterizes a large group of individual trajectories
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
4
Fokker-Plank Equation (KFE)
The Fokker-Plank equation gives
means for computing the probability
density p (t , x ) (mass of probabilities)
at any time.
Probability
for point xi
is useless
Probability
for range xi
has meaning
In this case of Wiener process the
Fokker-Plank equation can be solved
analytically as the Gaussian density
with zero mean and increasing in time
variance
1
p (t , x ) 
e
2 t

x2
2t
 N (0, t )( x )
Animate for your own pleasure:
http://en.wikipedia.org/wiki/Kolmogorov_forward_equation
This solution tells that the probability density is not defined for t  0 . Therefore, the
initial probability can’t be concentrated in a single point (as  ( x ) function does) but
in a wide area.
Otherwise, we should deal with distribution P (t , dx )  p (t , x )dx instead of density.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
5
Evolution of Probability Density
Fokker-Plank equation is solved numerically for a Wiener process
Start
p0 ( x )
Finish
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
As we expected in the
previous formula the
mean here don’t
change, only the
variance expands
progressively.
6
Probability Density for the different Initial Density
Also the initial density
is different now, but
this has a relatively
small effect on the
evolution process
compared to the
previous slide.
Start
p0 ( x )
Finish
The Fokker-Plank
equation works well
with ‘any’ initial
density.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
7
General Process: Fokker-Plank Equation (KFE)
If a stochastic process has both the drift and diffusion coefficients
dz  a (t , z (t ))dt  b(t , z (t ))dW
the Fokker-Plank equation reads
p 1  2 2


b (t , x) p (t , x)    a (t , x) p (t , x) 
2 
t 2 x
x
p (0, x)  p0 ( x)
Although this is a linear PDE, it is not a pleasant one computationally, because
the process coefficients (drift, diffusivity) are under the derivation operators and
that makes the solution process less stable numerically (we see this soon).
In the next Lecture, we consider a Backward Kolmogorov Equation (BKE) that
has a better structure.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
8
General Fokker-Plank Equation (KFE)
For multivariable (vector) process
z  Rd
dz  a(t , z (t ))dt  b(t , z (t ))dW
the Fokker-Plank equation reads
p 1 d 3  2
 
t 2 i , j 1 xi x j
bb (t, x)  
T
ij

p 
 ai (t, x ) p 
i 1 xi
d
p (0, x)  p0 ( x)
In practice, we solve this equation in 2D or 3D. Computation in higher
dimensions, say 4D, are impractical (too expensive computationally).
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
9
Fokker-Plank Equation in Matrix Notations – skip
The Fokker-Plank equation in matrix notations
p 1
   T  :  bb T (t , x) p    a (t , x) p 
t 2
is expressed through the Frobenius inner product*)
A:B 
d
T
T
A
B

tr
A
B

tr
AB
 ij ij
i , j 1
This is a component-wise multiplication of the matrix elements (thought as
they are vectors). The matrix B is a symmetric matrix:
bb (t , x)   bb (t , x)  ;
T
T
T

  vector;   
 x
T

y

 row

z 
*)
In linear algebra, the trace is sum of diagonal elements.
Here we have trace of a product that is similar to dot product of vectors.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
10
Boundary Condition
Theoretically, the problem should be solved in unbounded domain:   x   .
a  x  b.
However, in practice, it is solved in an essentially bounded domain:
Hence, we must insure that a solution p (t , x ) is fully inside
the domain using some heuristics. For example:
Computing tricks. At 1st we set the flux equal to zero on the boundary
p
x
xa
p

x
0
(Neumann condition – soft boundary)
x b
If we see that the solution is firmly inside the domain (and zero on the boundary),
then in the next step, we force solution to be zero on the boundary
p (t , a )  0
p (t , b)  0
(Dirichlet condition – severe boundary)
Frequently, the Dirichlet boundary gives a more stable result numerically than
the less firmly defined Neumann condition (more sensitive).
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
11
Probability Density is normalized – checking mean
If the initial distribution p0 ( x ) is a probability density then solution p (t , x ) of
KFE is also probability density. Integral over a probability density is one
(normalized). This gives a checking mean: compute the area under function
Q(t )   p(t , x)dx

and if this integral is one always,
Q(t )  1, then the solution is correct.
Before computing of moments we should insure that we have a probability density
p (t , x) and not just an unnormalized density c(t , x) .
In a ‘bad end’ we can ‘fix’ the problem by extra normalization – we replace an
unnormalized density with normalized (artificially) density
c(t , x)
c(t , x)  p (t , x) 
Q (t )
This is in what we sometimes run in moments computing.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
12
Computing of Mean and Variance
The Fokker-Plank equation can be solved numerically with Comsol solver and
the probability densities found. After having these densities, the interesting us
moments can be calculated by numerical integration that is simple to implement
under Comsol.
The mean Mz (t )  m(t ) and variance  (t ) of the process
dz  a (t , z (t ))dt  b(t , z (t ))dW
can be found by integration over the probability density

m(t ) 

xp (t , x)dx

 (t ) 


x 2 p(t , x )dx  m2 (t )

Next we apply these ideas to the geometric Brownian motion and then to the
portfolio problems.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
13
Characterization of Linear Process
In the case of geometric Brownian motion
dz = Az (t )dt  Bz (t )dW
the mean and variance are statistics of the log-normal distribution. They are not
the same good as in the Gaussian distribution statistics. Actually, we should
calculate more statistics (moments) that characterize asymmetry of the distribution
and so on. For calculation of the more moments we need the probability density
that can be found from solution of the KFE that is a simple equation in this case

2

 Bx  p(t , x )
1
p(t , x )

t
2
2
x 2
    Axp(t, x) 
x
If we expanded the KFE, it reads (spatially variable coefficients)
2
2
p (t , x)
p (t , x)
2 B  p (t , x )
2
2
x
2
B
A
B
x


 A  p (t , x)




2
x
2
t
x
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
14
Proof for training – skip

2

 Bx  p
1
2
2
x 2
B
2
   x p

2 x x
2
B 2   2 p

x
xp
2



2 x  x

p
p 
B2  2 2 p
2
2
2



x
x
p
x
2  x 2
x
x 
  Axp 
p
 Ax  Ap
x
x
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
15
Probability Density of Geometric Brownian motion
Start
p0 ( x )
Normalized
Finish
Process mean vale changes: 46
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
16
Probability Density for different Initial Density
Start
p0 ( x )
Finish
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
17
Stochastic Characterization
Q: What we did with the Fokker-Plank equation?
A: We mapped the initial distribution (e.g. initial wealth) in long-range future
distribution (now known as distribution of possible events).
6.5
6
5.5
5
4.5
4
3.5
3
0
0.2
0.4
0.6
Time, s
0.8
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
1
18
Characterization of Wealth (Investment Problem)
The Merton`s optimal investment model
dz   r+ opt (   r )  z (t )dt   opt z (t )dW



B
A
applied in the former Fokker-Plank equation, reads




2
p(t , x ) 1  


  x  p(t , x )     r+  opt (   r )  x p(t , x ) 
2  opt
 x  

2 x  
t

2
A


B


2
All the raw moments mn (t )  Mz n (t ) and central moments  n (t )  M  z (t )  m(t ) 
can be calculated by integration over the probability densities

mn (t ) 

n
x p (t , x)dx

Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
 n (t ) 

  x  m(t ) 
n
p (t , x) dx

19
n
Probability Density of Wealth – computing may fail reasons
Unfinished: 80% complete
Invested
p0 ( x )
Unfinished computing.
The computing process
was stopped beforehand
at 8 the terminal moment
at 10 is not reached.
Otherwise we will run into
the trouble on the left
boundary that is set too
narrow for FKE calculation
in the entire x-axis.
The presence of true
(financial) boundary
induces numerical
instability in the form of
vibration of the numerical
solution.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
It seemed this wasn’t a
good idea to stick to the
true (financial) boundary.
20
Probability Density of Wealth for different Initial Density
Unfinished: 80% complete
Invested
p0 ( x )
We don’t observe the same
numerical instability here.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
21
Characterization of Wealth (Investment-Consumption-Bequest)
The Merton’s optimal investment-consumption model
dz   r+ opt (   r )  z (t )dt  K (t ) z 2 (t )dt   opt z (t )dW



B
A
applied in the former Fokker-Plank equation, reads
p(t , x ) 1 

2
t
2
 
 x  p(t , x )
2
opt
x 2
    r+
opt
(   r )  K ( t ) x  xp ( t , x )
x
and if expanded, reads (temporarily and spatially variable coefficients)
2
2
p
p
2 B  p
2
2
x

x
2
B

A

K
(
t
)
x

B
 A  2 K (t ) x  p



2
2 x
t
x
For proof, see the next slide
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
22

Proof for training – skip
a (t , x ) = Ax  K (t ) x 2
b(t , x) = Bx
 2  b 2 (t , x) p 
x 2
2
2
    Bx  p    2 2 p

2
2 2  p
2 p
2
4
2
B
x
B
p
 
 2 B xp   B x


 x B
2


x 
x
x
x
x

 x 
  a (t , x ) p    Ax  K ( s ) x

x
x
2
p 
 Ax  K (s) x  px  p  A  2 K (s) x 
2
p B 2 x 2  2 p
p
2
2


x
2
B

A

K
(
s
)
x

B
 A  2 K ( s) x  p



2
t
2 x
x
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
23
Probability Density of Wealth with Bequest
Completed:10
Invested
p0 ( x )
In this case the
numerical instability
may be a problem.
Distribution
is unnormated
in this case (ocillation).
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
Kolmogorov backward
equation introduced later
is more stable
24
Probability Density of Wealth for different Initial Density
Completed
Invested
p0 ( x )
Also here we see some
numerical instability here.
It seemed, we must use the
backward Kolmogorov
equation instead
(Lecture 12)
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
25
Calculus of Moments
The numerical stability is not issue
always. With some other data the
probability density can be solved
smoothly as shown in figure above.
And then, the mean and standard
deviation can be computed are shown in
figure below.
Although a closed system for calculating
the statistical moments doesn’t exist, the
mean and variance can still be found by
integration over the probability densities.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
26
KFE Approximation with ODE: Representation Formula
The Chapman-Kolmogorov equation links together the probability density
p (t , x) with the transition density p (t , x; s, y ) and ‘initial’ density p ( s, y ) as the
sum over all possible events

p (t , x ) 
 p ( t , x ; s , y ) p ( s , y ) dy

p (t , x) – probability density for being in the state x at moment t
p ( s, y ) – ‘initial’ density (generalized for any instance before t )
p (t , x; s, y )  Pr  X (t )  x Y ( s )  y  – transition probability to overhand
a given state y at ‘initial’ moment s to the given state x at later moment t:
( s, y )  ( t , x )
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
27
Fokker-Plank Equation for Transition Probability
The transition probability p (t , x; s, y ) satisfies a modified Fokker-Plank equation
p(t , x; s, y ) 1  2 2


b ( s, y ) p ( t , x; s, y )    a ( s, y ) p ( t , x; s , y ) 
2 
t
2 x
x
lim p (t , x; s, y )   ( x)
t s
a ( s, y )  0 is simpler and again simpler in the case
of finite number of states xi : i  1... N , i.e. a Markov chain.
This equation without drift
Fokker-Plank equation for Markov chain is ODE
dP
 P(t )L
dt
P( s)  P( s)
P (0)  I
P – matrix of transition probabilities, L – matrix of events
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
28
Matrix P of Transition Probabilities
In Markov chain a process can jump from one state into other nearest state
with probability pij ( t )
Pr  X (t )  x j Y (t )  yi   pij (t )
in matrix representation this is
P (t )   pij (t ) : i, j  1,..., N 
Example: Population dynamics can be modeled as a Markov chain.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
29
Example: Population Dynamics
Four states xi  1, 2, 3, 4 represent the number of species per km2 that
changes suddenly at random moments.
This process depends on the birth
They define the matrix of events

and death

rates.
L
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
30
Matrix L of Events
 

L
0

0

0
0 
0 

  

    


0

 
Markov chain ( N   ) is simpler than corresponding continuous-time
process ( N   ). And also the ODE for Markov chain is simpler than the
corresponding Fokker-Plank equation (PDE)
dP
 P(t )L
dt
The matrix
P( s)  P( s)
P (0)  I
L of events is the remnant of the 2nd order derivative of KFE.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
31
Computing of Probabilities: Comsol 1D Settings
The Comsol diffusion-convection model (anisotropic and non-conservative)
c 
c
c
ts  ( D )  R  u
t x
x
x
1
can be put in the Fokker-Plank equation form
p 1  2 2


b (t , x) p    a (t , x) p 
2 
t 2 x
x
with following settings:
p  b 2 p
p

  bbx  a    bx2  bbxx  ax  p

 x 
2 x 
t x 
u
R
D
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
32
Proof for training – skip
 2  b2 p 
x 2
2

b
p    2 p


b 2   2 p   b 2 

 b
p
 p
 b

x x
x  x
x  x x x  x 
  b  2
 2 p p b 2
 2b 2  2 p
b p
 2b 

p 2  b
 2b
 2    b 2  p 
b
  x 

x x x x
x
x x
x x

x


 2 p
p
 2bbx
 2  bx2  bbxx  p
b
x x
x
  ap 
x
p
a
p
a  p
 a  ax p
x
x
x
p  b 2 p
p

  bbx  a    bx2  bbxx  ax  p

 x 
2 x 
t x 
u
R
D
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
33
Comsol Multivariable Settings
The Comsol diffusion-convection 3D model
p
ts  ( Dp )  R  up
t
1
can be put in the Fokker-Plank equation form (  is a vector)
p 1
   T  :  bb T (t , x) p    a (t , x) p 
t 2
with the following settings:
2D  bb T
R  p   T  : D + pa
u  ( T D ) T  a
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
34
Proof – skip
 ap   pa  a p
   :  Dp    Dp   p    : D  ( D)p
T
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
T
T
35
Comsol 2D Settings – skip (as far as you don’t need them)
2
 D
T
i , j 1
Dij


xi  x
   D11
y   D21
D12   D11 D21



D22   x
y
  2 D11
 x 2
2 2D
ij
T

D

:

 2
 

  D21
i , j 1 xi x j

 yx
2
2

b

b
2 D  bb T   11 12
b21b11  b22b12
 2 D12 
xy 
 2 D22 
2 
y 
b11b21  b12b22 

b212  b222 
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
D12 D22 

x
y 
 b11 b12 
b

b
b
 21 22 
36
Comsol Settings for KFE and moment computing
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
37
Orientation
In next Lecture we learn
Kolmogorov backward equation
This will be not a direct aiming gun but a more powerful gun. We face less
trouble in numerical solution in this case.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
38
Quiz Problem
Write the Fokker-Plank equation for computing of probability density in the case of
 Wiener process: dz = BdW
 geometric Brownian motion: dz = Az (t ) dt  Bz (t ) dW
 general process: dz  a ( t , z ( t ))dt  b( t , z ( t ))dW
How to convert the Fokker-Plank equation (above) in applicable Comsol model?
1
e
The probability density p (t , x ) is given (may be: p ( t , x ) 
2 t

x2
2t
). How to
compute unconditional average (mean) and variance? Just give the formulas.
In examination 1D case will be considered (only).
Half page for ‘Chiding Note’ is legal.
Prof. Robert Tenno, Control Engineering Laboratory, Helsinki University of Technology
39

Download Report