<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Processes
Introduction
In simulation and modelling we encounter a wide range of stochastic processes. But most fall into a few common categories
Step1:  ...and import some Python modules
Step2:  Ito processes
An Ito process is defined to be an adapted stochastic process that can be expressed as the sum of an integral with respect to a Wiener process and an integral with respect to time,
$$X_t = X_0 + \int_0^t \mu_s \, ds + \int_0^t \sigma_s \, dW_s,$$
or, in differential form,
$$dX_s = \mu_s \, ds + \sigma_s \, dW_s,$$
where $W$ is a Wiener process, $\sigma$ a predictable $W$-integrable process, $\mu$ predictable and Lebesgue-integrable. The integrability conditions can be expressed as
$$\int_0^t (\sigma_s^2 + |\mu_s|) \, ds < \infty.$$
$\mu$ and $\sigma$ are allowed to depend both on the time and current state, so we can write
$$X_t = X_0 + \int_0^t \mu(s, X_s) \, ds + \int_0^t \sigma(s, X_s) \, dW_s.$$
The function $\mu$ is referred to as drift, the function $\sigma$ as diffusion. The ItoProcess can thus be specified by providing these two functions
Step3:  It can then be approximated with a stochastic time discrete approximation, such as the Euler-Maruyama strong Taylor approximation scheme
Step4:  Since in this particular case the diffusion coefficient is constant, we could have defined X as
Step5:  Solved Ito processes
In case the stochastic differential equation
$$X_t = X_0 + \int_0^t \mu(s, X_s) \, ds + \int_0^t \sigma(s, X_s) \, dW_s$$
has a solution, it may be possible to compute the above integrals analytically. Although the rigorous theory of SDE requires one to specify what exactly one means by a "solution", we shall reserve the term solved Ito process for those Ito processes where these integrals can be computed analytically. Such processes should inherit from the SolvedItoProcess class and override its abstract method
Step6:  Given the time time0 and the process's value at that time, value0, and (if the process is stateful) the process's state, state0, at time0, as well as the random variate variate corresponding to the actual increment in the driving Brownian motion $W$, the propagate method will return the value of the process at time, time &gt;= time0. If propagate is implemented, there is no need to resort to approximate schemes, such as the Euler-Maruyama scheme demonstrated above.
In fact, the Ito process in our example above happens to be an Ornstein-Uhlenbeck process, whose solution is well known
Step7:  We make sure that we generate it with the same random seed...
Step8:  ...and verify that the graph is unchanged when we apply EulerMaruyama to this process, now instantiated as an OrnsteinUhlenbeckProcess, rather than an ItoProcess
Step9:  Instead of looping explicitly, we could have used the method run
Step10:  Now, since
Step11:  we don't need to apply Euler-Maruyama to produce a trajectory of this process and can use the propagate method instead
Step12:  Markov processes
Informally, a Markov process models the motion of a particle that moves around in a measurable space in a memoryless way. Such processes are extremely important and merit their own theory. The Wiener process and the Ornstein-Uhlenbeck process illustrated above are special cases. There also deserve special treatment in filtering theory &mdash; both Kalman and particle filtering.
In thalesians.tsa, Markov processes inherit from the abstract class MarkovProcess. Thus they provide the method
Step13:  This method represents the transition kernel of the Markov process
Step14:  Gaussian and Gauss-Markov processes
Gaussian processes play an important role in mathematical finance, machine learning, and in stochastic filtering, where the Kalman filter is the solution of this very special case of the filtering problem &mdash; the linear-Gaussian case.
Several Gauss-Markov processes are implemented in thalesians.tsa, notably the WienerProcess and OrnsteinUhlenbeckProcess, the latter being the only nontrivial stationary Gauss-Markov process. The Ornstein-Uhlenbeck process is ubiquitous in portfolio management and merits special treatment.
Specific processes
Wiener process
Univariate standard Wiener process
Step15:  Univariate variance-scaled Wiener process with drift
Step16:  Multivariate variance-scaled, correlated Wiener process with drift
Step17:  Brownian Bridge
Standard Brownian bridge
Step18:  Generalized Brownian bridge
Step19:  Multivariate Brownian bridge
Step20:  A more efficient method for simulating a Brownian bridge
Step21:  Geometric Brownian motion
Univariate geometric Brownian motion
Step22:  Multivariate variance-scaled, correlated geometric Brownian motion with drift
Step23:  Ornstein-Uhlenbeck process
Univariate Ornstein-Uhlenbeck process
Step24:  Multivariate Ornstein-Uhlenbeck process
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
Explanation: Processes
Introduction
In simulation and modelling we encounter a wide range of stochastic processes. But most fall into a few common categories: Ito processes, martingales, Markov processes, Gaussian processes, etc. We attempt to take this into account in our treatment of stochastic processes in thalesians.tsa, where we represent different categories of stochastic processes with distinct abstract data types.
Before we proceed, we need to enable Matplotlib to inline its graphs in this Jupyter notebook...
End of explanation
import os, sys
sys.path.append(os.path.abspath('../../main/python'))
import numpy as np
import matplotlib.pyplot as plt
import thalesians.tsa.numpyutils as npu
import thalesians.tsa.processes as proc
import thalesians.tsa.randomness as rnd
import thalesians.tsa.simulation as sim
Explanation: ...and import some Python modules:
End of explanation
X = proc.ItoProcess(drift=lambda t, x: -x, diffusion=lambda t, x: .25)
Explanation: Ito processes
An Ito process is defined to be an adapted stochastic process that can be expressed as the sum of an integral with respect to a Wiener process and an integral with respect to time,
$$X_t = X_0 + \int_0^t \mu_s \, ds + \int_0^t \sigma_s \, dW_s,$$
or, in differential form,
$$dX_s = \mu_s \, ds + \sigma_s \, dW_s,$$
where $W$ is a Wiener process, $\sigma$ a predictable $W$-integrable process, $\mu$ predictable and Lebesgue-integrable. The integrability conditions can be expressed as
$$\int_0^t (\sigma_s^2 + |\mu_s|) \, ds < \infty.$$
$\mu$ and $\sigma$ are allowed to depend both on the time and current state, so we can write
$$X_t = X_0 + \int_0^t \mu(s, X_s) \, ds + \int_0^t \sigma(s, X_s) \, dW_s.$$
The function $\mu$ is referred to as drift, the function $\sigma$ as diffusion. The ItoProcess can thus be specified by providing these two functions:
End of explanation
rnd.random_state(np.random.RandomState(seed=42), force=True);
ts = []; xs = []
for t, x in sim.EulerMaruyama(process=X, times=sim.xtimes(0., 100.)):
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: It can then be approximated with a stochastic time discrete approximation, such as the Euler-Maruyama strong Taylor approximation scheme:
End of explanation
X = proc.ItoProcess(drift=lambda t, x: -x, diffusion=.25)
Explanation: Since in this particular case the diffusion coefficient is constant, we could have defined X as
End of explanation
def propagate(self, time0, value0, time, variate=None, state0=None, random_state=None):
    raise NotImplementedError()
Explanation: Solved Ito processes
In case the stochastic differential equation
$$X_t = X_0 + \int_0^t \mu(s, X_s) \, ds + \int_0^t \sigma(s, X_s) \, dW_s$$
has a solution, it may be possible to compute the above integrals analytically. Although the rigorous theory of SDE requires one to specify what exactly one means by a "solution", we shall reserve the term solved Ito process for those Ito processes where these integrals can be computed analytically. Such processes should inherit from the SolvedItoProcess class and override its abstract method
End of explanation
X = proc.OrnsteinUhlenbeckProcess(transition=1, vol=.25)
Explanation: Given the time time0 and the process's value at that time, value0, and (if the process is stateful) the process's state, state0, at time0, as well as the random variate variate corresponding to the actual increment in the driving Brownian motion $W$, the propagate method will return the value of the process at time, time &gt;= time0. If propagate is implemented, there is no need to resort to approximate schemes, such as the Euler-Maruyama scheme demonstrated above.
In fact, the Ito process in our example above happens to be an Ornstein-Uhlenbeck process, whose solution is well known:
End of explanation
rnd.random_state(np.random.RandomState(seed=42), force=True);
Explanation: We make sure that we generate it with the same random seed...
End of explanation
ts = []; xs = []
for t, x in sim.EulerMaruyama(process=X, times=sim.xtimes(0., 100.)):
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: ...and verify that the graph is unchanged when we apply EulerMaruyama to this process, now instantiated as an OrnsteinUhlenbeckProcess, rather than an ItoProcess:
End of explanation
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, times=sim.xtimes(0., 100.))
df = sim.run(em)
plt.plot(df);
Explanation: Instead of looping explicitly, we could have used the method run:
End of explanation
isinstance(X, proc.SolvedItoProcess)
Explanation: Now, since
End of explanation
rnd.random_state(np.random.RandomState(seed=42), force=True);
x = 0.
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(1., 100., 1.), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: we don't need to apply Euler-Maruyama to produce a trajectory of this process and can use the propagate method instead:
End of explanation
def propagate_distr(self, time0, distr0, time):
    pass
Explanation: Markov processes
Informally, a Markov process models the motion of a particle that moves around in a measurable space in a memoryless way. Such processes are extremely important and merit their own theory. The Wiener process and the Ornstein-Uhlenbeck process illustrated above are special cases. There also deserve special treatment in filtering theory &mdash; both Kalman and particle filtering.
In thalesians.tsa, Markov processes inherit from the abstract class MarkovProcess. Thus they provide the method
End of explanation
def propagate(self, time0, value0, time, variate=None, state0=None, random_state=None):
    if self.noisedim != self.processdim:
        raise NotImplementedError('Cannot utilize the propagate_distr of the Markov process in propagate if noisedim != processdim; provide a custom implementation')
    if time == time0: return npu.tondim2(value0, ndim1tocol=True, copy=True)
    value0 = npu.tondim2(value0, ndim1tocol=True, copy=False)
    variate = npu.tondim2(variate, ndim1tocol=True, copy=False)
    distr = self.propagate_distr(time, time0, distrs.NormalDistr.creatediracdelta(value0))
    return distr.mean + np.dot(np.linalg.cholesky(distr.cov), variate)
Explanation: This method represents the transition kernel of the Markov process: given the marginal distribution at time0, distr0, and time &gt;= time0, the method returns the marginal distribution at time.
Solved Ito Markov processes
Processes that are both "solved Ito" (and therefore children of SolvedItoProcess) and Markov (and therefore children of MarkovProcess) should inherit from SolvedItoMarkovProcess. By default, their propagate is implemented in terms of their propagate_distr using the Dirac delta:
End of explanation
X = proc.WienerProcess()
x0 = 0.
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=1., step=1E-3))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = [x0]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 1., 1E-3), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Gaussian and Gauss-Markov processes
Gaussian processes play an important role in mathematical finance, machine learning, and in stochastic filtering, where the Kalman filter is the solution of this very special case of the filtering problem &mdash; the linear-Gaussian case.
Several Gauss-Markov processes are implemented in thalesians.tsa, notably the WienerProcess and OrnsteinUhlenbeckProcess, the latter being the only nontrivial stationary Gauss-Markov process. The Ornstein-Uhlenbeck process is ubiquitous in portfolio management and merits special treatment.
Specific processes
Wiener process
Univariate standard Wiener process
End of explanation
X = proc.WienerProcess(mean=3., vol=4.)
x0 = 7.
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=5., step=1E-3))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = [x0]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 5., 1E-3), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Univariate variance-scaled Wiener process with drift
End of explanation
X = proc.WienerProcess.create_from_cov(mean=[3., 5.], cov=[[16., -8.], [-8., 16.]])
x0 = npu.col(7., 8.)
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=5., step=1E-3))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = x0
ts = [0.]; xs = [x0.flatten()]
for t, v in zip(sim.xtimes(0., 5., 1E-3), rnd.multivariate_normals(ndim=2)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Multivariate variance-scaled, correlated Wiener process with drift
End of explanation
X = proc.BrownianBridge()
rnd.random_state(np.random.RandomState(seed=42), force=True);
ts = []; xs = []
for t, x in sim.EulerMaruyama(process=X, initial_value=0., times=sim.xtimes(0., 1., .005)):
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
rnd.random_state(np.random.RandomState(seed=42), force=True);
x = [0.]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 1., .005), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Brownian Bridge
Standard Brownian bridge
End of explanation
X = proc.BrownianBridge(10., 15., 0., 10.)
rnd.random_state(np.random.RandomState(seed=42), force=True);
ts = []; xs = []
for t, x in sim.EulerMaruyama(process=X, initial_value=10., times=sim.xtimes(0., 10., .005)):
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
rnd.random_state(np.random.RandomState(seed=42), force=True);
x = [10.]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 10., .005), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Generalized Brownian bridge
End of explanation
x0 = npu.col(10., 7.)
cov = [[1., -2.], [-2., 9.]]
X = proc.BrownianBridge.create_from_cov(x0, npu.col(15., 3.), 0., 10., cov)
rnd.random_state(np.random.RandomState(seed=42), force=True);
ts = []; xs = []
for t, x in sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(0., 10., .005)):
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
xs[-1]
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = x0
ts = [0.]; xs = [x0.flatten()]
for t, v in zip(sim.xtimes(0., 10., .005), rnd.multivariate_normals(ndim=2)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
xs[-1]
Explanation: Multivariate Brownian bridge
End of explanation
start_time = 0.
end_time = 10.
times = np.linspace(0., 10., 1000)
start_value = 10.
end_value = 15.
times_col = npu.to_ndim_2(times, ndim_1_to_col=True)
mean = start_value + (times_col - start_time) / (end_time - start_time) * (end_value - start_value)
cov = np.array([[(end_time - max(times[i], times[j])) * (min(times[i], times[j]) - start_time) / (end_time - start_time) for j in range(len(times))] for i in range(len(times))])
values = rnd.multivariate_normal(mean, cov)
all_times = np.concatenate(([start_time], times, [end_time]))
all_values = np.concatenate(([start_value], values, [end_value]))
plt.plot(all_times, all_values);
Explanation: A more efficient method for simulating a Brownian bridge
End of explanation
X = proc.GeometricBrownianMotion()
x0 = .3
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=1., step=1E-3))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = [x0]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 1., 1E-3), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
X = proc.WienerProcess()
x0 = .3
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=1., step=1E-3))
df = sim.run(em)
plt.plot(df);
Explanation: Geometric Brownian motion
Univariate geometric Brownian motion
End of explanation
X = proc.GeometricBrownianMotion.create_from_pct_cov(pct_drift=[3., 5.], pct_cov=[[16., -8.], [-8., 16.]])
x0 = npu.col(7., 8.)
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=1., step=1E-3))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = x0
ts = [0.]; xs = [x0.flatten()]
for t, v in zip(sim.xtimes(0., 1., 1E-3), rnd.multivariate_normals(ndim=2)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Multivariate variance-scaled, correlated geometric Brownian motion with drift
End of explanation
X = proc.OrnsteinUhlenbeckProcess(transition=1., vol=1.)
x0 = 0.
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=5., step=.01))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = [0.]
ts = [0.]; xs = [x]
for t, v in zip(sim.xtimes(0., 5., .01), rnd.multivariate_normals(ndim=1)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Ornstein-Uhlenbeck process
Univariate Ornstein-Uhlenbeck process
End of explanation
X = proc.OrnsteinUhlenbeckProcess.create_from_cov(
        transition=[[10., 0.], [0., 10.]],
        mean=[3., 5.],
        cov=[[9., -7.5], [-7.5, 25.]])
x0 = npu.col(7., 8.)
rnd.random_state(np.random.RandomState(seed=42), force=True)
em = sim.EulerMaruyama(process=X, initial_value=x0, times=sim.xtimes(start=0., stop=5., step=.01))
df = sim.run(em)
plt.plot(df);
rnd.random_state(np.random.RandomState(seed=42), force=True)
x = x0
ts = [0.]; xs = [x0.flatten()]
for t, v in zip(sim.xtimes(0., 5., .01), rnd.multivariate_normals(ndim=2)):
    x = X.propagate(ts[-1], x, t, v)
    ts.append(t); xs.append(x.flatten())
plt.plot(ts, xs);
Explanation: Multivariate Ornstein-Uhlenbeck process
End of explanation
<END_TASK>