Brownian.lhs

% The Existence of Brownian Motion
% Dominic Steinitz
% 24th September 2014

---
bibliography: Kalman.bib
---

In 1905, Einstein published five remarkable papers including, "Über
die von der molekularkinetischen Theorie der Wärme geforderte Bewegung
von in ruhenden Flüssigkeiten suspendierten Teilchen." which roughly
translates as "On the movement of small particles suspended in a
stationary liquid demanded by the molecular-kinetic theory of heat."
giving the first explanation of the phenomenon observed by Robert
Brown in 1827 of small particles of pollen moving at random when
suspended in water.

The eponymously named Brownian motion is a stochastic process
$\big(W_t\big)_{0 \le t \le 1}$ (presumably $W$ for Wiener) such that

1. $W_0(\omega) = 0$ for all $\omega$

2. For all $0 \le t_1 \le t \le t_2 \le, \ldots, t_n \le 1$, $W_{t_1},
W_{t_2} - W_{t_3}, \ldots W_{t_{n-1}} - W_{t_n}$ are independent.

3. $W_{t+h} - W_t \sim {\cal{N}}(0,h)$.

4. The map $t \mapsto W_t(\omega)$ is a continuous function of $t$ for
all $\omega$.

The
[Kolmogorov-Daniell](http://www.hss.caltech.edu/~kcb/Notes/Kolmogorov.pdf)
theorem guarantees that a stochastic process satisfying the first two
conditions exists but does not tell us that the paths are
continuous.

For example, suppose that we have constructed Brownian motion as
defined by the above conditions then take a random variable $U \sim
{\cal{U}}[0,1]$. Define a new process

$$
 \tilde{W}_t =
  \begin{cases}
   W & \text{if } t \neq U \\
   0 & \text{if } t = U
  \end{cases}
$$

This has the same finite dimensional distributions as $W_t$ as
$(W_{t_1}, W_{t_2}, W_{t_3}, \ldots W_{t_{n-1}}, W_{t_n})$ and
$(\tilde{W}_{t_1}, \tilde{W}_{t_2}, \tilde{W}_{t_3}, \ldots
\tilde{W}_{t_{n-1}}, \tilde{W}_{t_n})$ are equal unless $U \in \{t_1,
\ldots, t_n\}$ and this set has measure 0. This process satisifes
conditions 1--3 but is not continuous

$$
\lim_{t \uparrow U} \tilde{W}_t = \lim_{t \uparrow U} W_t = W_U
$$

and

$$
\mathbb{P}(W_U = 0) = \int_0^1
$$

Further this theorem is not constructive relying on the axiom of
choice. Most proofs of existence follow [@Ciesielski61]. Instead let
us follow [@liggett2010continuous].

It is tempting to think of Brownian motion as the limit in some sense
of a random walk.

$$
x_t^{(N)} = \sum_{i=1}^\floor{Nt}\frac{\xi_i}{N}
$$

where $0 \le t \le 1$ However, the processes, $x_t^{(1)}, x_t^{(2)},
\ldots$ are discontinuous and it is not clear how one would prove that
the limit of discontinuous processes is in fact continuous.

Let $\{\phi_i\}$ be a complete orthonormal basis for $L^2[0,1]$. That
is any $f$ for which $\int_0^1 f^2 \mathrm{d}\mu$ exists then

$$
f = \sum_{i=1}^\infty \langle f, \phi_i\rangle
$$

where $\mu$ is Lesbegue measure and $\langle\ldots,\ldots\rangle$ is
the inner product for $L^2$ defined as usual by

$$
\langle f, g\rangle = \int_0^1 fg\mathrm{d}\mu$
$$

We know such bases exist, for example, the [Fourier
expansion](http://en.wikipedia.org/wiki/Fourier_series) and [Legendre
polynomials](http://en.wikipedia.org/wiki/Legendre_polynomials).

We defined the so-called Haar wavelets for $n = 0, 1, \ldots$ and $k = 1, 3, 5, \ldots, 2^n - 1$.

$$
H_{n,k}(t) = +2^{(n-1)/2)} (k - 1)2^{-n} < t \le k2^{-n}
           = -2^{(n-1)/2)} k2^{-n}       < t \le (k + 1)2^{-n}
$$

Using Haskell's capabilities for dependently typed programming we can
express these as shown below.

> {-# OPTIONS_GHC -Wall                     #-}
> {-# OPTIONS_GHC -fno-warn-name-shadowing  #-}
> {-# OPTIONS_GHC -fno-warn-type-defaults   #-}
> {-# OPTIONS_GHC -fno-warn-unused-do-bind  #-}
> {-# OPTIONS_GHC -fno-warn-missing-methods #-}
> {-# OPTIONS_GHC -fno-warn-orphans         #-}

> {-# LANGUAGE DataKinds           #-}
> {-# LANGUAGE TypeOperators       #-}
> {-# LANGUAGE KindSignatures      #-}
> {-# LANGUAGE TypeFamilies        #-}
> {-# LANGUAGE TypeOperators       #-}
> {-# LANGUAGE Rank2Types          #-}
> {-# LANGUAGE ScopedTypeVariables #-}
> {-# LANGUAGE PolyKinds           #-}

> module Brownian where

> import GHC.TypeLits
> import Data.Proxy

> import Numeric.Integration.TanhSinh

> data Haar (a :: Nat) (b :: Nat) = Haar { unHaar :: Double -> Double }

> haar :: forall n k .
>         (KnownNat n, KnownNat k, (2 * k - 1 <=? 2^n - 1) ~ 'True) =>
>         Haar n (2 * k - 1)
> haar = Haar g where
>   g t | (k' - 1) * 2 ** (-n') < t && t <= k'       * 2 ** (-n') =  2 ** ((n' - 1) / 2)
>       | k'       * 2 ** (-n') < t && t <= (k' + 1) * 2 ** (-n') = -2 ** ((n' - 1) / 2)
>       | otherwise                                               =  0
>     where
>         n' = fromIntegral (natVal (Proxy :: Proxy n))
>         k' = 2 * (fromIntegral (natVal (Proxy :: Proxy k))) - 1

Now for example we can evaluate

> haar11 :: Double -> Double
> haar11 = unHaar (haar :: Haar 1 1)

    [ghci]
    haar11 0.75

but we if we try to evaluate *haar :: Haar 1 2* we get a type error.

> type family ZipWith (f :: a -> b -> c) (as :: [a]) (bs :: [b]) :: [c] where
>   ZipWith f (a ': as) (b ': bs) = (f a b) ': (ZipWith f as bs)
>   ZipWith f as        bs        = '[]

Rather than go too far into dependently typed programming which would
distract us from the existence proof, let us re-express this function
in a more traditional way (which will only give us errors at runtime).

> haarEtouffee :: Int -> Int -> Double -> Double
> haarEtouffee n k t
>   | n <= 0               = error "n must be >= 1"
>   | k `mod`2 == 0        = error "k must be odd"
>   | k < 0 || k > 2^n - 1 = error "k must be >=0 and <= 2^n -1"
>   | (k' - 1) * 2 ** (-n') < t && t <= k'       * 2 ** (-n') =  2 ** ((n' - 1) / 2)
>   | k'       * 2 ** (-n') < t && t <= (k' + 1) * 2 ** (-n') = -2 ** ((n' - 1) / 2)
>   | otherwise                                               =  0
>   where
>     k' = fromIntegral k
>     n' = fromIntegral n

Here are the first few Haar wavelets.

```{.dia height='300'}
import Brownian
import BrownianChart

dia = diag 2 "n = 1, k = 1" (xss!!0)
````

```{.dia height='300'}
import Brownian
import BrownianChart

dia = (diag 2 "n = 2, k = 1" (xss!!1) |||
       diag 2 "n = 2, k = 3" (xss!!2))
````

```{.dia height='300'}
import Brownian
import BrownianChart

dia = ((diag 2 "n = 3, k = 1" (xss!!3) |||
        diag 2 "n = 3, k = 3" (xss!!4) |||
        diag 2 "n = 3, k = 5" (xss!!5) |||
        diag 2 "n = 3, k = 7" (xss!!6)))
````

> schauderEtouffee :: Int -> Int -> Double -> Double
> schauderEtouffee n k t = result (absolute 1e-6 (parSimpson (haarEtouffee n k) 0 t))

> n :: Int
> n = 1000

> xss :: [[(Double, Double)]]
> xss = map (\(m, k) -> map (\i -> let x = fromIntegral i / fromIntegral n in (x, haarEtouffee m k x)) [0..n - 1]) [(1,1), (2,1), (2,3), (3,1), (3,3), (3,5), (3,7)]

> yss :: [[(Double, Double)]]
> yss = map (\(m, k) -> map (\i -> let x = fromIntegral i / fromIntegral n in (x, schauderEtouffee m k x)) [0..n - 1]) [(1,1), (2,1), (2,3), (3,1), (3,3), (3,5), (3,7)]

```{.dia height='300'}
import Brownian
import BrownianChart

dia = diag 0.5 "Foo" (yss!!0)
````

```{.dia height='300'}
import Brownian
import BrownianChart

dia = (diag 0.5 "Foo" (yss!!1) |||
       diag 0.5 "Baz" (yss!!2))
````

```{.dia height='300'}
import Brownian
import BrownianChart

dia = ((diag 0.5 "Foo" (yss!!3) |||
        diag 0.5 "Baz" (yss!!4) |||
        diag 0.5 "Baz" (yss!!5) |||
        diag 0.5 "Baz" (yss!!6)))
````

Bibliography and Resources
--------------------------