How to convert the parameters of multi-factors cheyette model (quasi-Gaussian model) from tenors to factors?

Question

The book "Interest Rate Modeling" by Andersen and Piterbarg is an extermely fascinating book on interest rate derivatives.
Recently, I have encoutered some issues while reading this book.
How to convert the parameters of multi-factors cheyette model from tenors to factors while keeping the variance of $S(t,T)$(Swap Rate) unchanged?
When calibrateing multi-factors cheyette model, we first calibrate the parameters for each tenor, similar to one-factor model.Then, we convert the parameters for multi-tenors into parameters for multi-factors. Right?
In the section 12.1.7 and 13.3.2 of the book "interest rate modeling", Andersen and Piterbarg tell us how to convert the parameters via restoring the variance/covariance of $f(t,T+ \delta_i)$. this is all well and good.
But my question is how to restore the variance of $S(t,T+\delta_i)$ ? Because the converted model must be able to Marked to Market.
In particular, the instataneous forward rate $$f(t,T) = f(0,T) + \left(x_t + yG(t,T)\right)\frac{h(T)}{h(t)} \tag{1}$$ $$x_t \sim N(0, \lambda^2\sigma^2 ) \tag{2}$$ let $D$ denotes the cholesky decomposed factors for correlation coefficent matrixs of $f(t,T+\delta_i)$. so we obtian $$\begin{pmatrix} \lambda_1\sigma_1&0&\cdots&0 \\\\ 0&\lambda_2\sigma_2&\cdots&0 \\\\ \vdots&\vdots&\ddots&\vdots \\\\ 0&0&\cdots&\lambda_n\sigma_n \\\\ \end{pmatrix}D= H_f \begin{pmatrix} \eta_{11}&\eta_{12}&\cdots&\eta_{1n} \\\\ \eta_{21}&\eta_{22}&\cdots&\eta_{2n} \\\\ \vdots&\vdots&\ddots&\vdots \\\\ \eta_{n1}&\eta_{n2}&\cdots&\eta_{nn} \\\\ \end{pmatrix}^T \tag{3}$$ where, the n×n matrix of $H^f$ denotes the factors's parameters. hence $$\begin{pmatrix} \eta_{11}&\eta_{12}&\cdots&\eta_{1n} \\\\ \eta_{21}&\eta_{22}&\cdots&\eta_{2n} \\\\ \vdots&\vdots&\ddots&\vdots \\\\ \eta_{n1}&\eta_{n2}&\cdots&\eta_{nn} \\\\ \end{pmatrix}^T = H^{-1}_f \begin{pmatrix} \lambda_1\sigma_1&0&\cdots&0 \\\\ 0&\lambda_2\sigma_2&\cdots&0 \\\\ \vdots&\vdots&\ddots&\vdots \\\\ 0&0&\cdots&\lambda_n\sigma_n \\\\ \end{pmatrix}D \tag{4}$$ So far, I have no questions. let's move on to simulate the $P(t,T)$. $$\ln P(t,T) =\ln P(0,T) -x_tG(t,T) - 0.5*G^T(t,T)yG(t,T) \tag{5}$$ where, $$G(t,T) = \int_t^T{\frac{h(s)}{h(t)}\text{d}s} \tag{6}$$ I found it is impossible to restore the variance of $P(t,T+\delta_i)$, not to mention restoring the variance of the Swap Rate. Is it ?

The following is my thought process for problem-solving:
The goal of the convertion is to restore the variance of $S(t,T+\delta_i)$, while keeping the correlation cofficient of $f(t,T+\delta_i)$ unchanged. And the correlation cofficient of $S(t,T+\delta_i)$ and the variance of $f(t,T+\delta_i)$ are unknown and not important.

Let's define the variance/covariance $$\Sigma_f =\begin{pmatrix} y_1&0&\cdots&0\\ 0&y_2&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&y_n\end{pmatrix} \rho \begin{pmatrix} y_1&0&\cdots&0\\ 0&y_2&\cdots&0\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&y_n\end{pmatrix} \tag{7}$$ Obviously, the correlation cofficient of $\Sigma_f$ is $\rho$.

Furthermore, accroding to formula(1) and (2), obtain $$\Sigma_f = H_f^T \Sigma_x H_f \tag{8}$$ Hence, $$\Sigma_x = (H_f^{-1})^T \Sigma_f H_f^{-1} = (H_f^{-1})^T Y\rho Y H_f^{-1}\tag{9}$$ where, $Y = diag\left(\begin{pmatrix} y_1&y_2&\cdots&y_n\end{pmatrix}^T \right)$.

And let denote the variance/covariance of $S(t, T+\delta_i)$ as $$\Sigma_s = M_s^T \Sigma_x M_s \tag{10}$$ where, $$M_s(t,T+\delta_i) = \frac{\partial S}{\partial x} (t,T+\delta_i) \tag{11}$$

Accroding to formula (9) and (10), obtain $$\Sigma_s = M_s^T (H_f^{-1})^T Y \rho Y H_f^{-1} M_s \tag{12}$$ Remember, we only need to restore the variance of Swap rate. Hence $$ diag\left(M_s^T (H_f^{-1})^T Y \rho Y H_f^{-1} M_s\right)= \begin{pmatrix} \sigma_{s,1}^2&\sigma_{s,2}^2& \cdots& \sigma_{s,n}^2\end{pmatrix}^T \tag{13}$$ let simplify the formula(13), obtain $$diag\left(A^T Z\rho Z A\right) = \begin{pmatrix} 1& 1&\cdots& 1\end{pmatrix}^T \tag{14}$$ Where,
$A= H_f^{-1}M_s$ is asymmetric invertible matrix,
$Z = diag\left(\begin{pmatrix} \frac{y_1}{\sigma_{s,1}}&\frac{y_2}{\sigma_{s,2}}&\cdots&\frac{y_n}{\sigma_{s,n}}\end{pmatrix}^T \right)$.

At last, my question is,
(a) how to solve the formula (14)? Does an analytic solution of $Z$ exist ?
(b) In order for $A^T Z\rho Z A$ be a correlation cofficient matrix, what conditions does $Z$ need to satisfy ?

Any help and hints are very much appreciated.

Answer 1

I've come up with a numerical solution method. Suppose $$A^TZ\rho ZA =\rho_S = \begin{pmatrix} 1&\cos(\theta_1 - \theta_2)&\cdots&\cos(\theta_1 - \theta_{n-1}) \\\\ \cos(\theta_2 - \theta_1)&1&\cdots&\cos(\theta_2 - \theta_{n-1}) \\\\ \vdots&\vdots&\ddots&\vdots \\\\ \cos(\theta_{n-1} - \theta_1)&\cos(\theta_{n-1} - \theta_2)&\cdots&1 \\\\ \end{pmatrix}_{n \times n} \tag{15}$$ where, $(\theta_1,...,\theta_{n-1})^T$ are the unknown parameters. Hence, $$\rho(\theta_1,...,\theta_{n-1})= (A^TZ)^{-1}\rho_S(ZA)^{-1}\tag{16}$$ Minimize the norm: $||\rho(\theta_1,...,\theta_{n-1})-\rho||$, to get result $(\theta_1,...,\theta_{n-1})$.
However, the solution obtained by the optimization generally has poor precision and is very time-consuming, which does not meet the requirements for computational accuracy and speed in exotic pricing. It is hoped that there can be a better solution.