ssvie 0.1
Numerical approximation for linear stochastic Volterra integral equations of the second kind
Requires Python: >=3.10

# ssvie

Numerical approximation for linear stochastic Volterra integral equations of the second kind, i.e.

$$X_t = f(t) + \int_0^t k_1(s,t) X_s ds + \int_0^t k_2(s,t) X_s dB_s \quad t \in [0,T),$$

where the function $f \colon [0,T) \to \mathbb{R}$ and the kernels $k_1, k_2 \colon [0,T)^2 \to \mathbb{R}$ are known, $X_t$ is an unknown process, $(B_t)_{t \in [0,T)}$ is a Brownian Motion and $\int_0^t k_2(s,t) X_s dB_s$ is the Itô-integral.

# 1. Block pulse functions (BPFs)

## 1.1 Defintion

Define the $m$-set of BPFs as

$$\phi_i(t) = \begin{cases} 1 & \text{, if } (i-1) \cdot h \leq t < i \cdot h \\ 0 & \text{, otherwise, }  \end{cases}$$

where $t \in [0,T)$, $i = 1, \ldots,m$ and $h = \frac{T}{m}$. The elementary properties of BPFs are

1. Disjointness:

   $$\phi_i(t) \phi_j(t) = \delta_{ij} \phi_i(t),$$

   where $\delta_{ij}$ is the Kronecker delta.
2. Orthogonality:

   $$\int_0^T \phi_i(t) \phi_j(t) dt = h \cdot\delta_{ij}$$
3. Completeness:

   $$\int_0^T f^2(t) dt = \sum_{i=0}^\infty f_i^2 \| \phi_i(t) \|^2,$$

   where

   $$f_i = \frac{1}{h} \int_0^T f(t) \phi_i(t) dt$$

   is the block pulse coefficient with respect to the $i$-th BPF $\phi_i(t)$.
4. Vector form:

   Let

   $$\phi(t) = ( \phi_1(t),\ldots, \phi_m(t))^T$$

   and

   $$F = (f_1,\ldots,f_m)^T.$$

## 1.2 Function approximation

A square integrable function $f \colon [0, T) \to \mathbb{R}$ can be expanded into a block pulse series

$$f(t) \simeq \hat{f}_m(t) = \sum_{i=1}^m f_i \phi_i(t),$$

where $f_i$ is the block pulse coefficient with respect to the $i$-th BPF $\phi_i(t)$. In vector form we have,

$$f(t) \simeq \hat{f}_m(t) = F^T \phi(t) = \phi(t)^TF,$$

where $F = (f_1,\ldots,f_m)^T$.

Let $k(s,t) \in L^2([0,T) \times [0,T))$. It can be similarly expanded with respect to BPFs such as

$$k(s,t) \simeq \hat{k}_m(s,t) = \phi^T(s) K \phi(t) = \phi^T(t) K \phi(s),$$

where $K=(k_{ij})_{i,j=1,\ldots,m}$ is the $m \times m$ block pulse coefficient matrix with

$$k_{ij} = \frac{1}{h^2} \int_0^T \int_0^T k(s,t) \phi_i(s)\phi_j(t) dt ds,$$

where $h=\frac{T}{m}$.

## 1.3 Integration operational matrix

It is

$$\int_0^t \phi(s) ds \simeq P \phi(t),$$

where

$$P = \frac{h}{2} \begin{pmatrix} 1 & 2 & 2 & \ldots & 2 \\
0 & 1 & 2 & \ldots & 2 \\ 
0 & 0 & 1 & \ldots & 2\\ 
\vdots & \vdots & \vdots & \ddots & \vdots \\ 
0 & 0 & 0 & \ldots & 1 \end{pmatrix}_{m \times m}$$

 is the so called operational matrix of integration.

# 2. Stochastic integration operational matrix

It is

$$\int_0^t f(s) dB_s \simeq \int_0^t F^T \phi(s) dB_s \simeq F^T P_S \phi(t),$$

where

$$P_S = \begin{pmatrix} B_{0.5h} & B_h & B_h & \ldots &  B_h\\
0 & B_{1.5h} - B_h & B_{2h} - B_h & \ldots & B_{2h} - B_h \\
0 & 0 & B_{2.5h} - B_{2h} & \ldots & B_{3h} - B_{2h}\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
0 & 0 & 0 & \ldots & B_{(m-0.5)h} - B_{(m-1)h} \end{pmatrix}_{m \times m}$$

is the so called stochastic operational matrix of integration.

# 3. Solving stochastic Volterra integral equations by using a stochastic operational matrix

Consider the linear stochastic Volterra integral equation of the second kind

$$X_t = f(t) + \int_0^t k_1(s,t) X_s ds + \int_0^t k_2(s,t) X_s dB_s \quad t \in [0,T),$$

where the function $f \colon [0,T) \to \mathbb{R}$ and the kernels $k_1, k_2 \colon [0,T)^2 \to \mathbb{R}$ are known, $X_t$ is an unknown process, $(B_t)_{t \in [0,T)}$ is a Brownian Motion and $\int_0^t k_2(s,t) X_s dB_s$ is the Itô-integral.

Then approximate $X_t$, $f(t)$, $k_1(s,t)$ and $k_2(s,t)$ as follows

+ $X_t \simeq X^T\phi(t) = \phi^T(t)X$
+ $f(t) \simeq F^T \phi(t) = \phi^T(t)F$
+ $k_1(s,t) \simeq \phi^T(s) K_1 \phi(t) = \phi^T(t) K_1^T \phi(s)$
+ $k_2(s,t) \simeq \phi^T(s) K_2 \phi(t) = \phi^T(t) K_2^T \phi(s)$


Set

$$B_1 = \frac{h}{2} \begin{pmatrix} k_11^{(1)} & 2k_{12}^{(1)} & 2k_{13}^{(1)} & \ldots & 2k_{1m}^{(1)} \\
0 & k_{22}^{(1)} & 2 k_{23}^{(1)} & \ldots & 2k_{2m}^{(1)} \\
0 & 0 & k_{33}^{(1)} & \ldots & 2k_{3m}^{(1)} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \ldots & k_{mm}^{(1)} \end{pmatrix}$$

and

$$B_2 = \begin{pmatrix} k_{11}^{(2)} B_{0.5h} & k_{12}^{(2)} B_h & k_{13}^{(2)} B_h & \ldots & k_{1m}^{(2)} B_h \\
0 & k_{22}^{(2)} (B_{1.5h} - B_h) & k_{23}^{(2)} (B_{2h} - B_h) & \ldots & k_{2m}^{(2)} (B_{2h} - B_h) \\
0 & 0 & k_{33}^{(2)} (B_{2.5h} - B_{2h}) & \ldots & k_{3m}^{(2)} (B_{3h} - B_{2h}) \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \ldots & k_{mm}^{(2)} (B_{(m-0.5)h} - B_{(m-1)h}) \end{pmatrix}.$$

A numerical solution to

$$X_t = f(t) + \int_0^t k_1(s,t) X_s ds + \int_0^t k_2(s,t) X_s dBS$$

is then given by the block pulse coefficient matrix $X$ of the unknown stochastic process $X_t$ which is obtained by solving

$$\left( I - B_1 - B_2 \right)^T X = F.$$

# 4. Examples

```{python}
import numpy as np
from ssvie import StochasticVolterraIntegralEquation
```

1. Consider

   $$X_t = 1 + \int_0^t s^2 X_s ds + \int_0^t sX_s dB_s \quad t \in [0,0.5)$$

   with the exact solution

   $$X_t = \exp \left( \frac{t^3}{6} + \int_0^t s dB_s \right)$$

   ```{python}
   # Define function and kernels
   def f_1(x):
     return 1

   def k1_1(s,t):
     return s**2

   def k1_2(s,t):
     return s

   example_1 = StochasticVolterraIntegralEquation(f_1, k1_1, k1_2, t=0.5)
   example_1_solution = example_1.solve_numerical(solve_method='bpf', m=100)
   ```
2. Consider

   $$X_t = \frac{1}{12} + \int_0^t \cos (s) X_s ds + \int_0^t \sin (s) X_s dB_s \quad t \in [0,0.5)$$

   with the exact solution


   $$X_t = \frac{1}{12} \exp \left( - \frac{t}{4} + \sin(t) + \frac{\sin (2t)}{8} + \int_0^t \sin (s) dB_s \right)$$

   ```{python}
   # Define function and kernels
   def f_2(x):
     return 1/12

   def k2_1(s,t):
     return np.cos(s)

   def k2_2(s,t):
     return np.sin(s)

   example_1 = StochasticVolterraIntegralEquation(f_2, k2_1, k2_2, t=0.5)
   example_1_solution = example_1.solve_numerical(solve_method='bpf', m=100)
   ```
3. Consider

   $$X_t = \frac{1}{3} + \int_0^t \ln (s+1) X_s ds + \int_0^t \sqrt{\ln (s+1)} X_s dB_s \quad t \in [0,0.5)$$

   with the exact solution

   $$X_t = \exp \left( - \frac{1}{2} t + \frac{1}{2} \ln (t+1) t + \frac{1}{2} \ln (t+1) + \int_0^t \sqrt{\ln (s+1)} dB_s \right)$$

   ```{python}
   def f_3(x):
     return 1

   def k3_1(s,t):
     return np.cos(s)

   def k3_2(s,t):
     return np.sin(s)

   example_3 = StochasticVolterraIntegralEquation(f_3, k3_1, k3_2, t=0.5)
   example_3_solution = example_3.solve_numerical(solve_method='bpf', m=100)
   ```

# Reference

(1) Deb, A., & Roychoudhury, S. (2018). "Control System Analysis and Identification with MATLAB®: Block Pulse and Related Orthogonal Functions". CRC Press.

(2) Maleknejad, K., Khodabin, M., & Rostami, M. (2012). "Numerical solution of stochastic Volterra integral equations by a stochastic operational matrix based on block pulse functions". Mathematical and computer Modelling, 55(3-4), 791-800. [https://doi.org/10.1016/j.mcm.2011.08.053](https://doi.org/10.1016/j.mcm.2011.08.053)