Time series generation for option pricing on quantum computers using tensor network

An option is originally a financial contract between two parties, one of which, the option holder, possesses the right, but not the obligation, to either buy or sell an underlying asset at a pre-determined price (strike), on a future date. If the holder has a right to buy (resp. sell) the underlying asset, the option is referred to as a call (resp. put) option. Let $S_{t}$ be an asset price at time $t$ and consider a call option written on that asset exercisable on a maturity date $t=T$. At the maturity date, the holder exercises a right to buy the asset if its price is higher than the strike $K$. This is equivalent to getting a profit $S_{T}-K$, since the asset is worth $S_{T}$. Conversely, if $S_{T}<K$, the holder chooses not to exercise the right and receives nothing. In total, this option is equivalent to the contract in which the option holder receives the payoff

$\displaystyle f_{\rm pay}(S_{T})=\max\{S_{T}-K,0\}\,,$

(1)

at $t=T$.

Derived from this kind of option called a European option, many kinds of options with various payoffs have been developed. In particular, our study focuses on path-dependent options, whose payoff depends not on the asset price at a single time point but on the values at multiple time points or the entire path of the asset price. In the following, we present examples of such options considered in this work.

The Heston model [35] is a popular mathematical model to describe the dynamics of the asset price and its volatility and has become widely used in the financial industry. It is one of stochastic volatility (SV) models, meaning it incorporates a stochastic process of the volatility as well as the asset price. SV models can be used when market prices of European options written on the asset do not match the prices given by the BS model, which is an often observed phenomenon called the volatility smile [41].

In the Heston model, the dynamics of the asset is given by the following SDE in the risk-neutral measure:

	$\displaystyle dS_{t}$	$\displaystyle=\sqrt{\nu_{t}}S_{t}dW_{t}^{S}$		(12)
	$\displaystyle d\nu_{t}$	$\displaystyle=\kappa(\theta-\nu_{t})dt+\xi\sqrt{\nu_{t}}dW_{t}^{\nu}$		(13)

where $W_{t}^{S}$ and $W_{t}^{\nu}$ are one-dimensional Wiener process with correlation $\rho$, which appears as $dW_{t}^{S}dW_{t}^{\nu}=\rho dt$, and $\kappa,\theta,\xi$ are positive constants. The Heston model describes the squared volatility $\nu_{t}$ as another stochastic process in Eq. (13), known as the CIR process [42], which exhibits mean-reverting. Note that the BS model corresponds to the case where $\nu_{t}$ is a positive constant, meaning the volatility is not dynamical, as opposed to the Heston model.

In this subsection, we briefly review option pricing with a quantum computer. The process of this algorithm is decomposed into three parts:

1. 1.

   Loading the probability distributions of price paths into a quantum state

2. 2.

   Encoding the payoff function of the option into the amplitude of an ancilla qubit

3. 3.

   Utilizing QAE to estimate the expectation value of the payoff function

As a first step, we prepare a quantum oracle $\mathcal{P}$, which encode the probability of price path $\{S_{t_{j}}\}_{j}=(S_{t_{0}},S_{t_{1}},\cdots S_{t_{M}})$ such that

$\displaystyle\mathcal{P}\ket{0}\rightarrow\sum_{\{S_{t_{j}}\}_{j}}\sqrt{p(\{S_{t_{j}}\}_{j})}\ket{\{S_{t_{j}}\}_{j}}\,,$

(14)

where $p(\{S_{t_{j}}\}_{j})$ is a probability associated with the path $\{S_{t_{j}}\}_{j}$ and $\ket{\{S_{t_{j}}\}_{j}}=\ket{S_{t_{0}}}\cdots\ket{S_{t_{M}}}$ is the computational basis state that represents the finite-precision binary representations of $S_{t_{0}},\cdots,S_{t_{M}}$. Then, in the second step, we operate the oracle $\mathcal{F}$, which acts on an ancilla qubit as

$\displaystyle\mathcal{F}\ket{\{S_{t_{j}}\}_{j}}\ket{0}\rightarrow\ket{\{S_{t_{j}}\}_{j}}\left(\sqrt{1-f_{\rm pay}(\{S_{t_{j}}\}_{j})}\ket{0}+\sqrt{f_{\rm pay}(\{S_{t_{j}}\}_{j})}\ket{1}\right)\,.$

(15)

As a result, the combination of oracles $\mathcal{Q}=\mathcal{F}\mathcal{P}$ creates a quantum state

$\displaystyle\sum_{\{S_{t_{j}}\}_{j}}\sqrt{p(\{S_{t_{j}}\}_{j})}\ket{\{S_{t_{j}}\}_{j}}\left(\sqrt{1-f_{\rm pay}(\{S_{t_{j}}\}_{j})}\ket{0}+\sqrt{f_{\rm pay}(\{S_{t_{j}}\}_{j})}\ket{1}\right)\,,$

(16)

which encodes the desired expectation value in the amplitude, as we see that the probability of measuring $\ket{1}$ in the last qubit is equal to

$\displaystyle\sum_{\{S_{t_{j}}\}_{j}}p(\{S_{t_{j}}\}_{j})f_{\rm pay}(\{S_{t_{j}}\}_{j})=E[f_{\rm pay}(\{S_{t}\})]\,.$

(17)

We do not cover the last step that uses QAE to elucidate the above probability; interested readers may see [8, 9]. By leveraging QAE, the query complexity to achieve the error $\epsilon$ scales as $O(1/\epsilon)$ while in classical algorithm it increases as $O(1/\epsilon^{2})$, which provides the quadratic speed-up of the Monte Carlo method. However, it is important to note that, as explained in Section 1, implementing the state-preparation operation $\mathcal{P}$ on quantum devices may be computationally demanding, and it is desired to find an efficient way to realize it.

In this study we propose employing MPS as a generative model $T_{\theta}$. To this end, we first need to discretize our dataset so that

$\displaystyle\mathcal{T}\rightarrow\bar{\mathcal{T}}=\{\bar{x}_{t}^{i}\}_{i=1}^{N},$

(18)

Here we transform continuous variables $x_{t}$ into integers in $\{0,\cdots,2^{m}-1\}$ corresponding to $m$-bit binary values as

$\displaystyle x_{j}\rightarrow\bar{x}_{j}=\left\lfloor(2^{m}-1)\times\frac{x_{j}-x_{\rm min}}{x_{\rm max}-x_{\rm min}}\right\rfloor$

(19)

with $x_{\rm min}:=\min_{i,j}x_{j}^{i}$ and $x_{\rm max}:=\max_{i,j}x_{j}^{i}$. This can naturally be fit in the computational basis of a quantum state. Then, the MPS $\Psi$ is introduced as a rank-$M$ tensor with each leg having dimension $2^{m}$, which associate $\bar{x}_{t}$ with a real number:

$\displaystyle\Psi(\bar{x}_{t})=\sum_{\alpha_{1}=0}^{D_{1}-1}\cdots\sum_{\alpha_{M}=0}^{D_{M}-1}A^{(1)}_{\bar{x}_{1}\alpha_{1}}A^{(2)}_{\alpha_{1}\bar{x}_{2}\alpha_{2}}\cdots A^{(M)}_{\alpha_{M-1}\bar{x}_{M}}\,,$

(20)

where $A^{(j)}_{\alpha_{j-1}\bar{x}_{j}\alpha_{j}}\in\mathbb{R}^{D_{i-1}\times 2^{m}\times D_{i}}$ is a rank-3 tensor, and $D_{i}\in\mathbb{N}_{>0}$ is referred to as bond dimension with $D_{0}=D_{M}=0$. The $j$-th entry $\bar{x}_{j}\in\{0,\cdots,2^{m}-1\}$ in $\bar{x}_{t}$ corresponds to the second leg of the $j$-th tensor, whose dimension $2^{m}$ are called physical dimension. Having the MPS $\Psi$, we calculate the probability distributions of $\bar{x}_{t}$ by

$\displaystyle p_{\theta}(\bar{x}_{t})=\frac{|\Psi(\bar{x}_{t})|^{2}}{Z},\,$

(21)

with the partition function $Z=\sum_{\bar{x}_{t}}\Psi(\bar{x}_{t})$.

It is worth noting that most of previous studies on generative modeling with MPS typically assign each digit of $\bar{x}$ into a different tensor, which means that, in our time series setting, each $\bar{x}_{j}$ is further decomposed into $\bar{x}_{j}=a_{j,0}+a_{j,1}\times 2+\cdots a_{j,m-1}\times 2^{m-1}$ with $a_{j,k}\in[0,1]$. Consequently, the model would provide us with the following structure

$\displaystyle\Psi(\bar{x}_{t})=\sum_{\tilde{\alpha}_{1},\cdots,\tilde{\alpha}_{Mm}}\tilde{A}^{(1)}_{a_{1,0}\tilde{\alpha}_{1}}\tilde{A}^{(2)}_{\tilde{\alpha}_{1}a_{1,1}\tilde{\alpha}_{2}}\cdots A^{(Mm)}_{\tilde{\alpha}_{Mm-1}a_{M,m-1}}\,.$

(22)

In such a formulation, the MPS model deals with correlation between $\bar{x}_{t-1}$ and $\bar{x}_{t}$ only though the first and last digits of them, which we suppose may result in poor expressive power for time series. This consideration motivates us to formulate our model as in Eq. (20).

In order to train our MPS model, we introduce the Kullback-leibler (KL) divergence as a measure to discern and evaluate the dissimilarity between probability distributions,

$\displaystyle D_{KL}\left(p_{\theta}(\bar{x}_{t})|\pi(\bar{x}_{t})\right)=\sum_{i=1}^{N}p_{\theta}(\bar{x}_{t}^{i})\ln\left(\frac{\pi(\bar{x}_{t}^{i})}{p_{\theta}(\bar{x}_{t}^{i})}\right)\,,$

(23)

where we denote the probability distribution of discretized random variable $\bar{x}_{t}$ as $\pi(\bar{x}_{t})$. Given that the KL divergence represents the distance between two distributions, our primary aim is to minimize Eq. (23). Minimizing the KL divergence can be achieved minimizing the following negative log-likelihood,

$\displaystyle\mathcal{L}=-\frac{1}{|\mathcal{T}|}\sum_{\bar{x}_{t}^{i}}\ln p_{\theta}(\bar{x}_{t}^{i})\,,$

(24)

which is equivalent to the KL divergence up to a constant. Thus, we adopt the negative log-likelihood as the objective function for training the MPS model.

The training of tensor components in the MPS is conducted though the gradient descent technique. While in usual manner of the gradient descent whole tensor elements are simultaneously updated, in this work we take advantage of the sweeping algorithm as in [29]. In this algorithm, each of adjacent tensor component pairs is iteratively updated while the others unchanged, and the bond dimension between the components is automatically adjusted, which is an advantage compared with the usual gradient descent approach. This is in fact similar to the density matrix renormalization group (DMRG) algorithm. In the following, we explain above sweeping algorithm in our setting.

Let us suppose that we are updating the $A^{(j)}$ and $A^{(j+1)}$. We first merge them into rank-4 tensors,

$\displaystyle A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}=\sum_{\alpha_{j}}A^{(j)}_{\alpha_{j-1}\bar{x}_{j}\alpha_{j}}A^{(j+1)}_{\alpha_{j}\bar{x}_{j+1}\alpha_{j+1}}\,.$

(25)

Then we compute the gradient of $\mathcal{L}$ with respect to the merged tensor $A^{(j,j+1)}$, which is given by

$\displaystyle\frac{\partial\mathcal{L}}{\partial A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}}=\frac{1}{|\mathcal{T}|Z}\frac{\partial Z}{\partial A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}}-\frac{2}{|\mathcal{T}|}\sum_{\bar{x}_{t}}\frac{1}{\Psi(\bar{x}_{t})}\frac{\partial\Psi(\bar{x}_{t})}{\partial A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}}.$

(26)

Using Eq. (26), we update the element of the merged tensor $A^{(j,j+1)}$ as

$\displaystyle A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}\rightarrow A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}-\eta\frac{\partial\mathcal{L}}{\partial A^{(j,j+1)}_{\alpha_{j-1}\bar{x}_{j}\bar{x}_{j+1}\alpha_{j+1}}}$

(27)

with a learning rate $\eta$. After updating the element of $A^{(j,j+1)}$, we perform the singular value decomposition for updated $A^{(j,j+1)}$:

$\displaystyle A^{(j,j+1)}=U\Lambda V^{\dagger}$

(28)

Here, $A^{(j,j+1)}$ is regarded as a $2^{m}D_{j-1}\times 2^{m}D_{j+1}$ matrix with the index pair $(\alpha_{j-1},\bar{x}_{j})$ specifying the row and $(\bar{x}_{j+1},\alpha_{j+1})$ specifying the column. $U$ and $V$ are orthogonal matrices with size $2^{m}D_{j-1}\times 2^{m}D_{j-1}$ and $2^{m}D_{j+1}\times 2^{m}D_{j+1}$, respectively. $\Lambda$ is a $2^{m}D_{j-1}\times 2^{m}D_{j+1}$ matrix in which the singular values of $A^{(j,j+1)}$ are diagonally arranged in descending order from the $(1,1)$th to $(D,D)$th entries, where $D=\min\{2^{m}D_{j-1},2^{m}D_{j+1}\}$, and the other entries are 0. We then truncate the singular values, that is, replace small entries in $\Lambda$ with 0. Supposing that $D_{j}^{\prime}$ singular values remain, we reach the approximation

$\displaystyle A^{(j,j+1)}$

$\displaystyle=\tilde{U}\tilde{\Lambda}\tilde{V}^{\dagger},$

(29)

where $\tilde{\Lambda}$ is the upper-left $D_{j}^{\prime}\times D_{j}^{\prime}$ block of $\Lambda$, and $\tilde{U}\in\mathbb{R}^{2^{m}D_{j-1}\times D_{j}^{\prime}}$ (resp. $\tilde{V}\in\mathbb{R}^{2^{m}D_{j+1}\times D_{j}^{\prime}}$) consists of the first $D_{j}^{\prime}$ columns of $U$ (resp. $V$). Finally, we update the original rank-3 matrices by

$\displaystyle A^{(j)}=\tilde{U}\sqrt{\tilde{\Lambda}}\,,\quad A^{(j+1)}=\sqrt{\tilde{\Lambda}}\tilde{V}^{\dagger},$

(30)

where $A^{(j)}\in\mathbb{R}^{D_{j-1}\times 2^{m}\times D_{j}^{\prime}}$ is regarded as a $2^{m}D_{j-1}\times D_{j}^{\prime}$ matrix with $(\alpha_{j-1},\bar{x}_{j})$ specifies the row and $\alpha_{j}$ specifies the column, and $A^{(j+1)}\in\mathbb{R}^{D_{j}^{\prime}\times 2^{m}\times D_{j+1}^{\prime}}$ is regarded as a $D_{j}^{\prime}\times D_{j+1}2^{m}$ matrix with $\alpha_{j}$ specifies the row and $(x_{j+1},\alpha_{j+1})$ specifies the column. Note that the bond dimension $D_{j}$ has changed to $D_{j}^{\prime}$.

We start the algorithm with performing this procedure for the rightmost pair. After that, we move on to the next left pair, and repeat the same. When reaching the leftmost tensor, we reverse the updating process from left to right. The whole sweeping loop starting and ending at the rightmost tensor defines one epoch of the training. We then run some epochs to get the final result. The bond dimensions $D_{j}$ are adjusted as explained above, with some criterion for singular value truncation and an upper bound $D_{\rm max}$. Consult [29] for further details.

While our final goal is to prepare a quantum state encoding a probability distribution, one can generate samples from the MPS model. Thanks to the fact that the partition function of MPS can be calculated efficiently, it enables us to draw samples directly without help of such as Gibbs sampling. To begin with, we take the leftmost tensor in MPS and compute the following marginal probability:

$\displaystyle P(\bar{x}_{1})=\frac{\sum_{\bar{x}_{2},\cdots,\bar{x}_{M}=0}^{2^{m}-1}|\Psi(\bar{x}_{1},\bar{x}_{2},\cdots,\bar{x}_{M})|^{2}}{Z}\,,$

(31)

where we only contract the second and below physical legs of MPS in the numerators. According to $P(\bar{x}_{1})$, we can draw a sample of $\bar{x}_{1}$. Given that sample, the drawn probability of the next one $\bar{x}_{2}$ is given as the conditional probability:

$\displaystyle P(\bar{x}_{2}|\bar{x}_{1})=\frac{P(\bar{x}_{2},\bar{x}_{1})}{P(\bar{x}_{1})}\,,$

(32)

where $P(\bar{x}_{2},\bar{x}_{1})$ can similarly be calculated as $P(\bar{x}_{1})$. Using this conditional probability $P(\bar{x}_{2}|\bar{x}_{1})$, $\bar{x}_{2}$ can be drawn. Repeating this procedure one-by-one, we are able to obtain the series of samples ${\bar{x}_{t}}=(\bar{x}_{1},\bar{x}_{2},\cdots,\bar{x}_{n})$. Then, each integer $\bar{x}_{j}$ is converted to the real value $x_{j}$ by

$\displaystyle x_{j}=x_{\rm min}+\frac{\bar{x}_{j}}{2^{m}-1}(x_{\rm max}-x_{\rm min}),$

(33)

which can be regarded as an approximation of the original random variable with rounding error of order $O(2^{-m})$.

We here summarize the setting of our numerical experiments. We set these parameters as $\kappa=1.0$, $\theta=0.04$, $\xi=2$, and $\rho=-0.7$. The initial values of the asset price $S_{0}$ and the squared volatility $v_{0}$ are set to 100 and 0.04, respectively. To generate paths in the Heston model, we adopt the simple Euler-Maruyama discretization scheme [43], that is, we discretize the original stochastic process as below:

	$\displaystyle S_{i+1}$	$\displaystyle=S_{i}+\sqrt{\nu_{i}}S_{i}\sqrt{\Delta t}Z^{S}\,,$
	$\displaystyle\nu_{i+1}$	$\displaystyle=V_{i}+\kappa(\theta-\nu_{i})\Delta t+\xi\sqrt{\nu_{i}}\sqrt{\Delta t}Z^{\nu}\,,$		(34)

where $\Delta t=t_{i+1}-t_{i}$ is a time step, and $Z^{S}$ and $Z^{\nu}$ are random variables drown from the bivariate standard normal distribution with correlation $\rho$. To avoid a negative variance, we assume the reflection positivity, $\nu_{i}=-\nu_{i}$ if $\nu_{i}<0$. In the experiment, we set the time step as $\Delta t=1/250$, corresponding daily observations, and the length of time series as $M=5$. With these settings, we then generate $N=10000$ paths of the Heston model.

In the training of the MPS model, the physical dimension are varied as $2^{m}$ with $m=4,5,6$. That means we discretize each price into $m$-bit binary values. Concretely, we adopt the aforementioned way,

$\displaystyle S_{t}\rightarrow\bar{S_{t}}=\left\lfloor(2^{m}-1)\times\frac{S_{t}-S_{\rm min}}{S_{\rm max}-S_{\rm min}}\right\rfloor$

(35)

with $S_{\rm min}$ and $S_{\rm min}$ are the minimum and maximum of the asset price, respectively, over sample paths and $n$ time points. Since our training procedure allows us to optimize bond dimensions, we fix the maximum value of bond dimensions $D_{\rm max}$ in the training. In our numerical experiment, we consider different values of $D_{\rm max}=64,100,150$.

After the training, we generate $N=10000$ paths using the resultant model, converting the model output integers to real values as Eq. (33). Employing the Monte Carlo method, we utilize generated price paths to price a European call option and path-dependent options described in 2.1, Asian, lookback, up-and-out barrier call options. The payoffs in these options, which are originally defined with the entire path in continuous time, are now replaced with the discrete-time versions:

$\displaystyle f_{\rm pay}(S_{t_{1}},\cdots,S_{t_{M}})=\max\left\{\frac{1}{M}\sum_{j=1}^{M}S_{t_{j}}-K,0\right\}$

(36)

for an Asian option,

$\displaystyle f_{\rm pay}(S_{t_{1}},\cdots,S_{t_{M}})=\max\left\{\max_{j=1,\cdots,M}S_{t_{j}}-K,0\right\}$

(37)

for a lookback option, and

$\displaystyle f_{\rm pay}(S_{t_{1}},\cdots,S_{t_{M}})=\begin{cases}\max\left\{S_{t_{M}}-K,0\right\}&(\max_{j=1,\cdots,M}S_{t_{j}}\,<B)\\ 0&(\mathrm{otherwise})\end{cases}.$

(38)

Strikes of these options are set to $K=100$. As for the up-and-out Barrier option, we set the barrier level $B=105$. We also compute prices of these options using the paths in the original Heston model generated by Eq. (34) as a benchmark. By comparing the results in the two ways, we can analyze the effectiveness and accuracy of the MPS model in pricing these path-dependent options.

To compare results of the Heston model and the MPS model more closely, we also introduce the implied volatility (IV) $\sigma_{BS}(K,T)$ as follows:

$\displaystyle C=C_{BS}(S_{0},K,T,\sigma_{BS}(K,T))$

(39)

where $C$ is a price of a European call option with strike $K$ and maturity $T$. That is, given the market price $C$ we reversely calculate $\sigma_{BS}(K,T)$ from Eq. (39). With $C$ being the market price of the option, we can regard the IV as the market’s expectation of the future volatility of the underlying asset. One reason why the IV is widely used in practice is that it is a strike-independent indicator of the option price level: the option prices for different strikes largely differs but the IVs are comparable. It is also common to compare different option pricing methods in terms of IV.

We first consider the effect of physical dimensions in MPS models. For that purpose, we fix $D_{\rm max}=150$ and vary the physical dimensions $2^{m}$ with $m=4,5,6$. Figure 1 shows the resulted probability distributions of $S_{t_{1}}$, $S_{t_{3}}$, and $S_{t_{5}}$ for each physical dimension. There, we plot empirical distributions of $S_{t}$ for $N=10000$ samples generated by the Heston and MPS models. These figures indicate that the higher $m$ becomes the closer the distribution is to the Heston model. Note that, for $m=4$, which means prices are discretized into only $2^{4}=16$ grid points, the resulting distribution is localized around several values, as in the left graph in Figure 1, but such a phenomenon is mitigated for larger $m$.

We then evaluate the price of options, and the IV calculated from the value of the European call option, showing the result in Table 1. Similarly to the probability distributions at each time step, with higher physical dimensions, option prices as well as the IV calculated by the MPS model get closer to those by the Heston model. On the European option, the IV by the MPS model with $m=6$ has a difference of several tenths of a percent from that by the Heston model, which we can say is practically a good fit for a pricing model of path-dependent options. Also for the Asian and lookback options, the prices by the MPS model are close to those by the Heston model with differences similar to that of the European option. Thus, we can conclude that the MPS model with sufficiently large physical dimension successfully learns the Heston model, at least for the purpose of pricing these options. On the other hand, for the barrier option, the difference between the prices of the two models is still larger even at $m=6$. We may improve the MPS model by, e.g., taking larger $m$ or using tensor network architectures other than MPS, which can be considered in future works.

The MPS model has the other hyperparameter, which would be another factor for the model’s capability, that is, the maximal value of the bond dimension $D_{\rm max}$. In order to evaluate the effect of $D_{\rm max}$, we thus consider the different $D_{\rm max}$ with fixed $m$. In the below, we set $m=6$. Figure 2 shows the probability distributions of $S_{t}$ for various $D_{\rm max}$. It is observed that with higher $D_{\rm max}$ the distributions are closer to those of the Heston model. We also investigate the effect of $D_{\rm max}$ in terms of option prices and the IV as shown in Table 2, where the results of the Heston and MPS models with $D_{\rm max}=150$ are same as in Table 1. We find that, basically, the prices by the MPS model get closer to those by the Heston model as we increase $D_{\rm max}$. Only for the barrier option, the price difference becomes larger as $D_{\rm max}$ increases, which implies that the MPS model underestimates its value possibly because the MPS model generates more paths that excess the barrier level, which largely affects the conditional expectation value for the barrier option, than the Heston model. The outliers in generated price paths can cause the discrepancy as well. These problems may be resolved by adopting some regularization technique to the model during training and generation, or by using tensor network structures other than MPS, which should be addressed in the future work.