The purpose is to simulate a functional autoregressive process of the form
$$
X_t(u)=\sum_{k=1}^p \int_0^1\Psi_k(u,v) X_{t-k}(v)dv+\varepsilon_t(u),\quad 1\leq t\leq n.
$$
Here we assume that the observations lie in a finite dimensional subspace of the function space spanned by
Fourier basis functions \(\boldsymbol{b}^\prime(u)=(b_1(u),\ldots, b_d(u))\). That is \(X_t(u)=\boldsymbol{b}^\prime(u)\boldsymbol{X}_t\), \(\varepsilon_t(u)=\boldsymbol{b}^\prime(u)\boldsymbol{\varepsilon}_t\) and \(\Psi_k(u,v)=\boldsymbol{b}^\prime(u)\boldsymbol{\Psi}_k \boldsymbol{b}(v)\). Then it follows that
$$
\boldsymbol{X}_t=\boldsymbol{\Psi}_1\boldsymbol{X}_{t-1}+\cdots+ \boldsymbol{\Psi}_p\boldsymbol{X}_{t-p}+\boldsymbol{\varepsilon}_t.
$$
Hence the dynamic of the functional time series is described by a VAR(\(p\)) process.
In this mathematical model the law of \(\boldsymbol{\varepsilon}_t\) is determined by noise.
The matrices Psi[,,k] correspond to \(\boldsymbol{\Psi}_k\). If op.norms is provided, then the coefficient
matrices will be rescaled, such that the Hilbert-Schmidt norms of \(\boldsymbol{\Psi}_k\) correspond to the vector.