from aff_poly_sig.exp_sig import expsig, withwords, expsig_withwords, moments
Affine and polynomial methods for signatures
Documentation
The code documentation is available at https://sarasvaluto.github.io/AffPolySig/.
The original theoretical results can be found in the following papers.
C. Cuchiero, G. Gazzani, J. Möller, and S. Svaluto-Ferro. Joint calibration to SPX and VIX options with signature-based models. ArXiv e-prints, 2023. https://arxiv.org/abs/2301.13235.
C. Cuchiero, S. Svaluto-Ferro, and J. Teichmann. Signature SDEs from an affine and polynomial perspectiv. ArXiv e-prints, 2023. http://arxiv.org/abs/2302.01362
How to install
After cloning the repository, cd
(change directory) to the repo and enter this into your terminal:
pip install -e .
How to use
Expected signature and moments of a polynomial process
The goal of this code is to compute the expected signature \(\mathbb E[\mathbb X_T]\) of a polynomial process \(X\). In the one dimensional case it also provides an expression for its moments. As a first step import the following functions.
Next, define the parameters of the polynomial process of interest. Given the drift vector \[b(X_t)^i=b_i+\sum_{j=0}^db_{ij}X_t^j,\] and the diffusion matrix \[a(X_t)^{ij}=a_{ij}+\sum_{k=0}^da_{ijk}X_t^k+\sum_{k,h=0}^da_{ijkh}X_t^kX_t^h,\] we use the following parametrisation of the characteristics \[\begin{align*} b_{const}[i]&=b_i,& b_{lin}[i,j]&=b_{ij},\\ a_{const}[i,j]&=a_{ij},& a_{lin}[i,j,k]&=a_{ijk},& a_{quad}[i,j,k,h]=a_{ijkh}. \end{align*}\] Coefficients need then to be saved in a tuple as illustrated in the code below. Define also the initial condition \(x_0\).
For this example we consider a one dimensional Jacobi process without drift setting \(b(X_t)^0=0\) and \(a(X_t)^{00}=X_t^0(1-X_t)^0\). We set \(x_0=1/2\).
import numpy as np
#Dimension of the process
=1
dim
#Coefficients of the characteristics
=np.zeros(dim)
b_const=np.zeros((dim,dim))
b_lin=np.zeros((dim,dim))
a_const=np.zeros((dim,dim,dim))
a_lin=np.zeros((dim,dim,dim,dim))
a_quad
0]=1
a_lin[0]=-1
a_quad[
=(b_const,b_lin,a_const,a_lin,a_quad)
coeff
=np.zeros(dim)
x00]=1/2 x0[
The last parameters to define are given by the lenght of the expected signature we would like to compute (len_max
) and the time at which we would like to do it (T
).
=10
len_max=1 T
To get the expected signature we can then use the function expsig
.
expsig(coeff,x0,len_max,dim,T)
array([1.00000000e+00, 0.00000000e+00, 7.90150699e-02, 0.00000000e+00,
1.45583443e-03, 0.00000000e+00, 1.14548260e-05, 0.00000000e+00,
4.94622744e-08, 0.00000000e+00, 1.34438690e-10])
To obtain a better readable output one can use the function withwords
.
=expsig(coeff,x0,len_max,dim,T)
E withwords(E,dim)
[[1.0, []],
[0.07901506985356974, [0, 0]],
[0.0014558344297645905, [0, 0, 0, 0]],
[1.145482602482784e-05, [0, 0, 0, 0, 0, 0]],
[4.94622743536099e-08, [0, 0, 0, 0, 0, 0, 0, 0]],
[1.3443868989326117e-10, [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]]
The same result can be obtained directly using expsig_withwords
, which combines the two steps above.
expsig_withwords(coeff,x0,len_max,dim,T)
[[1.0, []],
[0.07901506985356974, [0, 0]],
[0.0014558344297645905, [0, 0, 0, 0]],
[1.145482602482784e-05, [0, 0, 0, 0, 0, 0]],
[4.94622743536099e-08, [0, 0, 0, 0, 0, 0, 0, 0]],
[1.3443868989326117e-10, [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]]
In one dimension (dim
=1) we can also use the function moments
to compute the moments of the process.
moments(coeff,x0,len_max,dim,T)
array([1. , 0.5 , 0.40803014, 0.36204521, 0.33448524,
0.31613774, 0.30305083, 0.29324856, 0.28563369, 0.27954838,
0.2745743 ])
Laplace transform in the Brownian setting
The goal of this code is to compute the Laplace transform \(\mathbb E[e^{\langle \mathbf u^{sig}, \mathbb X_t^{sig}\rangle}]\) and \(\mathbb E[e^{\langle \mathbf u^{pow}, \mathbb X_t^{pow}\rangle}]\), where \(X\) denotes a Brownian motion and \(\mathbb X^{sig}\) and \(\mathbb X^{pow}\) the corresponding extensions. Precisely, we consider the following two extensions: the signature \[\mathbb X_t^{sig}:=(1,X_t,\frac {X_t^2}2,\ldots),\] and the power sequence \[\mathbb X_t^{pow}:=(1,X_t,X_t^2,\ldots).\] As a first step import the following functions.
from aff_poly_sig.riccati_bm import appr_exp_sig, appr_exp_pow, MC, CoD
Next, we introduce the parameters of interest. We in particular have T
for the time horizon, K_u
for the lenght of (the approximation of) \(\mathbf u^{sig}\) and \(\mathbf u^{pow}\), u_sig
for \(\mathbf u^{sig}\) with the signature’s extension, and u_pow
for \(\mathbf u^{pow}\) with the power sequence’s extension. As an example we compute \[\mathbb E[\exp(-2\frac{X_t^2}2)]=\mathbb E[\exp(-X_t^2)]\] for each \(t\in[0,T]\), where \(X\) denotes a Brownian Motion and \(T=1\). The corresponding parameters are given in the following cell.
import math
import matplotlib.pyplot as plt
=1 # [0,T] time horizon
T
#u in terms of the sig lift
=5
K_u=np.zeros(K_u)
u_sig2]=-2
u_sig[
#u in terms of the powers lift
=u_sig.copy()
u_powfor k in range(0,K_u):
=u_pow[k]/math.factorial(k) u_pow[k]
One then just needs to fix the computational parameters, namely the grid for the computation of the signature (timegrid
) and the trucation’s level K
for the solution.
=1000
n_time= np.linspace(0,T,n_time)
timegrid =30 K
The desired Laplace transform can then be computed using appr_exp_sig
and appr_exp_pow
. The finals sig
and pow
denote the employed extensions.
=appr_exp_sig(u_sig,timegrid,K)
L_sig=appr_exp_pow(u_pow,timegrid,K) L_pow
To test the results, one can compare them to the results obtained by a Monte Carlo approximation. That’s what the MC
function does.
=100000 # number of samples
N=1000 # number of times ticks
n_MC= MC(u_sig,T,n_MC,N) MonteCarlo
#Monte Carlo
#MonteCarlo_CoD = CoD(np.real(MonteCarlo),n_MC, n_time)
0,T,n_MC),MonteCarlo,'c',label='MC');
plt.plot(np.linspace(
#Sig Lift
'm--',label='Sig');
plt.plot(timegrid,L_sig,
#Pow Lift
'g:',label='Pow');
plt.plot(timegrid,L_pow,
min(MonteCarlo.real-0.01),max(MonteCarlo.real)+0.01)
plt.ylim("t")
plt.xlabel(f'Solving the Riccati for u_sig={u_sig.real} with an ODE solver')
plt.title(; plt.legend()