This notebook explains our design decisions for simulating stochastic dynamical systems.


```python
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from kalman import LDSNet
```

# Weiner Process

Describes a stochastic process $W$ which is driven by noise $w$ and governed by

$$dW/dt=v$$

The process is characterized $W(t+\Delta t)-W(t)\sim\mathcal{N}(0, \Delta t)$ (i.e., normally-distributed with mean 0 and variance $\Delta t$.

Discretize the stochastic differential equation SDE by chopping up time into $\Delta t$-sized bins and apply the Euler step

\begin{align}
\frac{W_{t+\Delta t}-W_t}{\Delta t} &= v_t \\
W_{t+\Delta t} &= W_t + \Delta tv_t \\
\end{align}

So that $W_{t+\Delta t}-W_t~\sim\mathcal{N}(0, \Delta t)$, we'll need to set $v_t~\sim\mathcal{N}\left(0, \frac{1}{\Delta t}\right)$

As an alternative, we could use the update $W_{t+\Delta t} = W_t + \sqrt{\Delta t}\tilde{v}_t$ where $\tilde{v}_t\sim\mathcal{N}(0, 1)$



```python
T = 1
Q = 1 # Target variance per unit time
dts = [1, 0.5, 0.1, 0.05, 0.01][::-1]
runs = 100000
data = {}
for dt in dts:
    N = int(np.round(T/dt))
    dW = np.random.randn(runs, N)*np.sqrt(dt)*np.sqrt(Q)
    W = np.zeros((runs, N+1))
    W[:, 1:] = np.cumsum(dW, axis=1)
    data[dt] = W
```

Showing a sample paths for each dt


```python
# plot trajectories
N_PLOT = 2
for dt in dts:
    N = int(np.round(T/dt))
    trange = np.arange(N+1)*dt
    color = plt.plot(trange, data[dt][0, :], label=dt)[0].get_color()
    for plt_idx in range(N_PLOT)[1:]:
        plt.plot(trange, data[dt][plt_idx+1, :], color=color)
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")
plt.xlabel("Time")
plt.ylabel("W")
```

Note how regardless of time discretization, the variance grows at the same rate 


```python
# plot variance over time
plt.figure()
for dt in dts:
    N = int(np.round(T/dt))
    trange = np.arange(N+1)*dt
    var = np.var(data[dt], axis=0)
    plt.plot(trange, var, label=dt)
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")
plt.xlabel("Time")
plt.ylabel("Variance")
    
plt.figure()
for dt in dts:
    plt.hist(data[dt][:, -1], bins=100, density=True, histtype="step", label=str(dt))
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")
plt.title("Distribution of W at T=1")
plt.xlabel("W(1)")
```

## A Note of Caution

If we naively used $W_{t+\Delta t} = W_t + \Delta t\tilde{v}_t$ where $\tilde{v}_t\sim\mathcal{N}(0, 1)$, then the rate of variance growth over time is linear with $\Delta t$.


```python
for dt in dts:
    N = int(np.round(T/dt))
    dW = np.random.randn(runs, N)*dt
    W = np.zeros((runs, N+1))
    W[:, 1:] = np.cumsum(dW, axis=1)
    data[dt] = W

# plot trajectories
plt.figure()
N_PLOT = 2
for dt in dts:
    N = int(np.round(T/dt))
    trange = np.arange(N+1)*dt
    color = plt.plot(trange, data[dt][0, :], label=dt)[0].get_color()
    for plt_idx in range(N_PLOT)[1:]:
        plt.plot(trange, data[dt][plt_idx+1, :], color=color)
plt.xlabel("Time")
plt.ylabel("W")
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")

# plot variance over time
plt.figure()
for dt in dts:
    N = int(np.round(T/dt))
    trange = np.arange(N+1)*dt
    var = np.var(data[dt], axis=0)
    plt.plot(trange, var, label=dt)
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")
plt.xlabel("Time")
plt.ylabel("Variance")
    
plt.figure()
for dt in dts:
    plt.hist(data[dt][:, -1], bins=100, density=True, histtype="step", label=str(dt))
plt.legend(loc="center left", bbox_to_anchor=(1, 0.5), title="dt")
plt.title("Distribution of W at T=1")
plt.xlabel("W(1)")
```


```python

```

ead5991ce2bfc2a4e87064f3775cef46038b8813
ipynb
Jupyter Notebook
stochastic_differential_equation.ipynb
samfok/kalman
f26c190f110a4b0044b2d6dd9acb72a32dbe4532
MIT
2019-08-01T09:28:23.000Z
2019-08-01T09:28:23.000Z
stochastic_differential_equation.ipynb
samfok/kalman
f26c190f110a4b0044b2d6dd9acb72a32dbe4532
MIT
stochastic_differential_equation.ipynb
samfok/kalman
f26c190f110a4b0044b2d6dd9acb72a32dbe4532
MIT
2019-08-01T09:28:22.000Z
2020-11-20T09:48:31.000Z
Qwen/Qwen-72B
1. YES
2. YES