**2. Model Description**

Next, we define the model in a filtered probability space (<sup>Ω</sup>, F,P, F) where F0 contains all subsets of the (P−) null sets of F and F is right-continuous. We first provide the processes under the historical measure P, then followed by the processes under a (conveniently chosen) risk-neutral measure Q. Suppose that *X*(*t*)=( *<sup>X</sup>*1(*t*), ... , *Xn*(*t*)) is a vector of asset prices with the following P-measure representation:

$$\begin{split} \frac{dX\_i(t)}{X\_i(t)} &= \left\{ L\_i + c\_i \sum\_{j=1}^p a\_{ij}^2 \left( \sqrt{v\_j(t)} + \frac{b\_j}{\sqrt{v\_j(t)}} \right)^2 - \sum\_{j=1}^n \beta\_{ij} \ln(X\_j(t)) + \tilde{c}\_i \left( \sqrt{\tilde{v}\_i(t)} + \frac{\tilde{b}\_i}{\sqrt{v\_i(t)}} \right)^2 \right\} dt \\ &+ \sum\_{j=1}^p a\_{ij} \left( \sqrt{v\_j(t)} + \frac{b\_j}{\sqrt{v\_j(t)}} \right) d\mathcal{W}\_j(t) + \left( \sqrt{\tilde{v}\_i(t)} + \frac{\tilde{b}\_i}{\sqrt{\tilde{v}\_i(t)}} \right) d\tilde{\mathcal{W}}\_i(t) \\ d\upsilon\_j(t) &= \quad a\_j(\theta\_j - v\_j(t)) dt + \xi\_j \sqrt{v\_j(t)} dB\_j^P(t), \; j = 1, \ldots, p \\ d\bar{\upsilon}\_i(t) &= \quad \bar{\mathfrak{x}}\_i(\theta\_i - \bar{\upsilon}\_i(t)) dt + \xi\_i \sqrt{\bar{\upsilon}\_i(t)} dB\_i^P(t), \; i = 1, \ldots, n \end{split}$$

<sup>1</sup> This can also be applied to volatility indexes, such as those reported by the Chicago Board Options Exchange (CBOE), which are clearly a mean-reverting asset class with stochastic volatility.

The quadratic variation structure is \$*dBPj* (*t*), *dWj*(*t*)%= *<sup>ρ</sup>jdt*, \$*dB*#*Pi* (*t*), *dW*#*i*(*t*)%= *ρ*#*idt* and zero otherwise. In the language of factor analysis, *aij* is the *ij*th entry of the matrix of factor loadings (*A*) that captures the correlations among assets. The communalities are represented by *Vj*(*t*) = *vj*(*t*) + √*bj vj*(*t*) 2 (in matrix form, <sup>Λ</sup>*nxp* = *A* diag *V*1/2(*t*)) and the intrinsic residual variance is Ψ = diag(*V* # (*t*)), with *V* # *j*(*t*) = *v*˜*i*(*t*) + ˜ √*bi v*˜*i*(*t*) 2. This leads to a factors decomposition of the quadratic variation of asset prices as follows:

$$
\Sigma(t)dt = \left(\Lambda\Lambda' + \Psi\right)dt = \left(A\operatorname{diag}(V(t))A' + \operatorname{diag}(\tilde{V}(t))\right)dt
$$

Whenever necessary, we assume *n* = *p* and *A* = (*aij*)*n*×*p* to be an orthogonal matrix. In this setting, *ci* and *c*˜*i* represent risk premiums of asset *Xi*(*t*) associated with the common and intrinsic factors, respectively. *β* = (*βij*)*n*×*n* is an invertible matrix, which captures the spillover at the expected return level *Xi*(*t*) on asset *Xj*(*t*). In other words, it represents the impact from other assets on the long term average price of the current one.

Based on the quadratic variation relationship defined in this model, if we assume that *BPj* , *BPj* (*t*)<sup>⊥</sup>, *B* ˜ *P i* (*t*), *B*˜ *Pi* (*t*)<sup>⊥</sup> are independent Brownian motions with −1 ≤ *ρj* ≤ 1 and −1 ≤ *ρ*˜*i* ≤ 1. Then,

$$d\mathcal{W}\_{\boldsymbol{j}}(t) = \rho\_{\boldsymbol{j}} d\boldsymbol{B}\_{\boldsymbol{j}}^{P}(t) + \sqrt{1 - \rho\_{\boldsymbol{j}}^{2}} d\boldsymbol{B}\_{\boldsymbol{j}}^{P}(t)^{\perp}$$

$$d\widetilde{\mathcal{W}}\_{\boldsymbol{i}}(t) = \rho\_{\boldsymbol{i}} d\boldsymbol{B}\_{\boldsymbol{i}}^{P}(t) + \sqrt{1 - \rho\_{\boldsymbol{i}}^{2}} d\boldsymbol{B}\_{\boldsymbol{i}}^{P}(t)^{\perp}.$$

*vj j* = 1, ... , *n* and *v*˜*i i* = 1, ... , *n* follow standard CIR processes, hence *<sup>α</sup>j*, *θj*, and *ξj* are positive constants satisfying *<sup>α</sup>jθj* ≥ *ξ*2*j*2 (the Feller condition). Similarly, *<sup>α</sup>*˜*i*, ˜*θi*, and ˜*ξi* are positive constants satisfying *α*˜*i* ˜ *θi* ≥ ˜ *ξi*2 2 . Note that the Feller condition in CIR model guarantees that the process remains positive.

The transformation *Y* = ln *X* would create a multivariate Ornstein–Uhlenbeck process with a 4/2 stochastic factor structure:

$$\begin{split} d\mathcal{Y}\_{l}(t) &= \left\{ \mathcal{L}\_{l} + \left(c\_{l} - \frac{1}{2}\right) \sum\_{j=1}^{n} a\_{lj}^{2} \left(\sqrt{\upsilon\_{j}(t)} + \frac{b\_{j}}{\sqrt{\upsilon\_{j}(t)}}\right)^{2} - \sum\_{j=1}^{n} \mathcal{S}\_{l} \mathcal{Y}\_{j}(t) + \left(\varepsilon\_{l} - \frac{1}{2}\right) \left(\sqrt{\mathcal{J}\_{l}(t)} + \frac{\mathcal{S}\_{l}}{\sqrt{\upsilon\_{l}(t)}}\right)^{2} \right\} dt \\ &+ \sum\_{j=1}^{n} a\_{lj} \left(\sqrt{\upsilon\_{j}(t)} + \frac{b\_{j}}{\sqrt{\upsilon\_{j}(t)}}\right) dW\_{l}(t) + \left(\sqrt{\mathcal{J}\_{l}(t)} + \frac{\mathcal{S}\_{l}}{\sqrt{\mathcal{J}\_{l}(t)}}\right) d\bar{W}\_{l}(t) \end{split}$$

This general model includes a notable particular case, which is a direct generalization of Grasselli (2017) to a factor setting when *βij* = 0, *i*, *j* = 1, ... , *n* ; this case is studied in more detail and it is named **FG**, given its analytical flexibility.
