Consistent Re-Calibration of the Discrete-Time Multifactor Vasi\v{c}ek Model

The discrete-time multifactor Vasi\v{c}ek model is a tractable Gaussian spot rate model. Typically, two- or three-factor versions allow one to capture the dependence structure between yields with different times to maturity in an appropriate way. In practice, re-calibration of the model to the prevailing market conditions leads to model parameters that change over time. Therefore, the model parameters should be understood as being time-dependent or even stochastic. Following the consistent re-calibration (CRC) approach, we construct models as concatenations of yield curve increments of Hull-White extended multifactor Vasi\v{c}ek models with different parameters. The CRC approach provides attractive tractable models that preserve the no-arbitrage premise. As a numerical example, we fit Swiss interest rates using CRC multifactor Vasi\v{c}ek models.


Introduction
The tractability of affine models, such as the Vasiček [2] and the Cox-Ingersoll-Ross [3] models, has made them appealing for term structure modeling. Affine term structure models are based on a (multidimensional) factor process, which in turn describes the evolution of the spot rate and the bank account processes. No-arbitrage arguments then provide the corresponding zero-coupon bond prices, yield curves and forward rates. Prices in these models are calculated under an equivalent martingale measure for known static model parameters. However, model parameters typically vary over time as financial market conditions change. They may, for instance, be of a regime switching nature and need to be permanently re-calibrated to the actual financial market conditions. In practice, this re-calibration is done on a regular basis (as new information becomes available). This implies that model parameters are not static and, henceforth, may also be understood as stochastic processes. The re-calibration should preserve the no-arbitrage condition, which provides side constraints in the re-calibration. The aim of this work is to discuss these side constraints with the help of the discrete-time multifactor Vasiček interest rate model, which is a tractable, but also flexible model. We show that re-calibration under the side constraints naturally leads to Heath-Jarrow-Morton [4] models with stochastic parameters, which we call consistent re-calibration (CRC) models [1]. These models are attractive in financial applications for several reasons. In risk management and in the current regulatory framework [5], one needs realistic and tractable models of portfolio returns. Our approach provides tractable non-Gaussian models for multi-period returns on bond portfolios. Moreover, stress tests for risk management purposes can be implemented efficiently in our framework by selecting suitable models for the parameter process. While an in-depth market study of the performance of CRC models remains to be done, we provide in this paper some evidence of improved fits. The paper is organized as follows. In Section 2, we introduce Hull-White extended discrete-time multifactor Vasiček models, which are the building blocks for CRC in this work. We define CRC of the Hull-White extended multifactor Vasiček model in Section 3. Section 4 specifies the market price of risk assumptions used to model the factor process under the real-world probability measure and the equivalent martingale measure, respectively. In Section 5, we deal with parameter estimation from market data. In Section 6, we fit the model to Swiss interest rate data, and in Section 7, we conclude. All proofs are presented in Appendix A. Choose a (sufficiently rich) filtered probability space (Ω, F, F, P * ) with discrete-time filtration F = (F(t)) t∈N0 , where t ∈ N 0 refers to time point t∆. Assume that P * denotes an equivalent martingale measure for a (strictly positive) bank account numeraire (B(t)) t∈N0 . B(t) denotes the value at time t∆ of an investment of one unit of currency at Time 0 into the bank account (i.e., the risk-free rollover relative to ∆). We use the following notation. Subscript indices refer to elements of vectors and matrices. Argument indices refer to time points. We denote the n × n identity matrix by 1 ∈ R n×n . We also introduce the vectors 1 = (1, . . . , 1) ∈ R n and e 1 = (1, 0, . . . , 0) ∈ R n .

Discrete-Time Multifactor Vasiček Model
We choose n ∈ N fixed and introduce the n-dimensional F-adapted factor process: X = (X(t)) t∈N0 = (X 1 (t), . . . , X n (t)) t∈N0 , which generates the spot rate and bank account processes as follows: where t ∈ N 0 ; empty sums are set equal to zero. The factor process X is assumed to evolve under P * according to: with initial factor X(0) ∈ R n , b ∈ R n , β ∈ R n×n , Σ 1 2 ∈ R n×n and (ε * (t)) t∈N = (ε * 1 (t), . . . , ε * n (t)) t∈N being F-adapted. The following assumptions are in place throughout the paper. Assumption 1. We assume that the spectrum of matrix β is a subset of (−1, 1) n and that matrix Σ 1 2 is non-singular. Moreover, for each t ∈ N, we assume that ε * (t) is independent of F(t − 1) under P * and has standard normal distribution ε * (t) P * ∼ N (0, 1).
Remark. In Assumption 1, the condition on matrix β ensures that 1 − β is invertible and that the geometric series generated by β converges. The condition on Σ 1 2 ensures that Σ = Σ 1 2 (Σ 1 2 ) is symmetric positive definite. Under Assumption 1, Equation (2) defines a stationary process; see [6], Section 11.3. (1) and (2) is called the discrete-time multifactor Vasiček model. Under the above model assumptions, we have for m > t:

The model defined by Equations
Remark. For m > t, the conditional distribution of X(m), given F(t), depends only on the value X(t) at time t∆ and on lag m−t. In other words, the factor process (2) is a time-homogeneous Markov process.
At time t∆, the price of the zero-coupon bond (ZCB) with maturity date m∆ > t∆ with respect to filtration F and equivalent martingale measure P * is given by: For the proof of the following result, see Appendix A.
Theorem 2. The ZCB prices in the discrete-time multifactor Vasiček Models (1) and (2) with respect to filtration F and equivalent martingale measure P * have an affine term structure: with A(m − 1, m) = 0, B(m − 1, m) = 1∆ and for m − 1 > t ≥ 0: In the discrete-time multifactor Vasiček Models (1) and (2), the term structure of interest rates (yield curve) takes the following form at time t∆ for maturity dates m∆ > t∆: with the spot rate at time t∆ given by Y (t, t + 1) = 1 X(t) = r(t).

Hull-White Extended Discrete-Time Multifactor Vasiček Model
The possible shapes of the Vasiček yield curve (4) are restricted by the choice of the parameters b ∈ R n , β ∈ R n×n and Σ ∈ R n×n . These parameters are not sufficiently flexible to exactly calibrate the model to an arbitrary observed initial yield curve. Therefore, we consider the Hull-White extended version (see [7]) of the discrete-time multifactor Vasiček model. We replace the factor process defined in (2) as follows. For fixed k ∈ N 0 , let X (k) satisfy: with starting factor X (k) (k) ∈ R n , e 1 = (1, 0, . . . , 0) ∈ R n and function θ : N → R. Model assumption (5) corresponds to (2), where the first component of b is replaced by the time-dependent coefficient (b 1 + θ(i)) i∈N and all other terms ceteris paribus. Without loss of generality, we choose the first component for this replacement. Note that parameter b 1 is redundant in this model specification, but for didactical reasons, it is used below. The time-dependent coefficient θ is called the Hull-White extension, and it is used to calibrate the model to a given yield curve at a given time point k∆. The upper index (k) denotes that time point and corresponds to the time shift we apply to the Hull-White extension θ in Model (5). The factor process X (k) generates the spot rate process and the bank account process as in (1). The model defined by (1,5) is called the Hull-White extended discrete-time multifactor Vasiček model. Under these model assumptions, we have for m > t ≥ k: Remark. For m > t ≥ k, the conditional distribution of X (k) (m), given F(t), depends only on the factor X (k) (t) at time t∆. In this case, factor process (5) is a time-inhomogeneous Markov process. Note that the upper index (k) in the notation is important since the conditional distribution depends explicitly on the lag m − k.
Theorem 3. The ZCB prices in the Hull-White extended discrete-time multifactor Vasiček model (1,5) with respect to filtration F and equivalent martingale measure P * have affine term structure: In the Hull-White extended discrete-time multifactor Vasiček model (1,5), the yield curve takes the following form at time t∆ for maturity dates m∆ > t∆ ≥ k∆: with spot rate at time t∆ given by Y (k) (t, t + 1) = 1 X (k) (t).
Remark. Note that the coefficient B(t, m) in Theorem 3 is not affected by the Hull-White extension θ and depends solely on m − t, whereas the coefficient A (k) (t, m) depends explicitly on the Hull-White extension θ.

Calibration of the Hull-White Extended Model
We consider the term structure model defined by the Hull-White extended factor process X (k) and calibrate the Hull-White extension θ ∈ R N to a given yield curve at time point k∆. We explicitly introduce the time index k in Model (5) because the CRC algorithm is a concatenation of multiple Hull-White extended models, which are calibrated at different time points k∆, see Section 3 below. Assume that there is a fixed final time to maturity date M ∆ and that we observe at time k∆ the yield curve y(k) ∈ R M for maturity dates (k + 1)∆, . . . , (k + M )∆. For these maturity dates, the Hull-White extended discrete-time multifactor Vasiček yield curve at time k∆, given by Theorem 3, reads as: For given starting factor X (k) (k) ∈ R n and parameters b ∈ R n , β ∈ R n×n and Σ ∈ R n×n , our aim is to choose the Hull-White extension θ ∈ R N such that we get an exact fit at time k∆ to the yield curve y(k), that is, The following theorem provides an equivalent condition to (7), which allows one to calculate the Hull-White extension θ ∈ R N explicitly.
Theorem 4 shows that the Hull-White extension can be calculated by inverting the (M − 1) × (M − 1) lower triangular positive definite matrix C(β).

Consistent Re-Calibration
The crucial extension now is the following: we let parameters b, β and Σ vary over time, and we recalibrate the Hull-White extension in a consistent way at each time point, that is according to the actual choice of the parameter values using Theorem 4. Below, we show that this naturally leads to a Heath-Jarrow-Morton [4] (HJM) approach to term structure modeling.

Consistent Re-Calibration Algorithm
Assume that (b(k)) k∈N0 , (β(k)) k∈N0 and (Σ(k)) k∈N0 are F-adapted parameter processes with β(k) and Σ(k) satisfying Assumption 1, P * -a.s., for all k ∈ N 0 . Based on these parameter processes, we define the n-dimensional F-adapted CRC factor process X , which evolves according to Steps (i)-(iv) of the CRC algorithm described below. Thus, factor process X will define a spot rate model similar to (1).
In the CRC algorithm, Steps 3.1.1-3.1.3 below are executed iteratively.

3.1.2.
Increments of the Factor Process from k → k+1. Assume factor X (k), parameters b(k), β(k) and Σ(k) and Hull-White extension θ (k) are given. Define the Hull-White extended model X (k) = (X (k) (t)) t≥k by: with starting value X (k) (k) = X (k), F(k)-measurable parameters b(k), β(k) and Σ(k) and Hull-White extension θ (k) . We update the factor process X at time (k + 1)∆ according to the X (k) -dynamics, that is, we set: This provides F(k + 1)-measurable yield curve at time (k + 1)∆ for maturity dates m∆ > (k + 1)∆: , and recursively for m − 1 > t ≥ k: This is exactly the no-arbitrage price under P * if the parameters b(k), β(k) and Σ(k) and the Hull-White extension θ (k) remain constant for all t > k.
(End of algorithm.) Remark. For the implementation of the above algorithm, we need to consider the following issue. Assume we start the algorithm at Time 0 with initial yield curve y(0) ∈ R M . At times k∆, for k > 0, calibration of θ (k) ∈ R M −1 requires yields with times to maturity beyond M ∆. Either yields for these times to maturity are observable, and the length of θ (k) is reduced in every step of the CRC algorithm or an appropriate extrapolation method beyond the latest available maturity date is applied in every step.

Heath-Jarrow-Morton Representation
We analyze the yield curve dynamics (Y(k, ·)) k∈N0 obtained by the CRC algorithm of Section 3.1. Due to re-calibration (10), the yield curve fulfills the following identity for m > k + 1: where the first line is based on the F(k)-measurable parameters (b(k), β(k), Σ(k)) and Hull-White extension θ (k) , and the second line is based on the F(k + 1)-measurable parameters and Hull-White extension (b(k + 1), β(k + 1), Σ(k + 1), θ (k+1) ) after CRC Step (iii). Note that in the re-calibration only (b(k + 1), β(k + 1), Σ(k + 1)) can be chosen exogenously, and the Hull-White extension θ (k+1) is used for consistency property (10). Our aim is to express Y(k + 1, m) as a function of X (k) and Y(k, m). Using Equations (9) and (11), we have for m > k + 1: This provides the following theorem; see Appendix A for the proof.
Theorem 5. Under equivalent martingale measure P * , the yield curve dynamics (Y(k, ·)) k∈N0 obtained by the CRC algorithm of Section 3.1 has the following HJM representation for m > k + 1: Key observation. Observe that in Theorem 5, a remarkable simplification happens. Simulating the CRC algorithm (9) and (10) to future time points k∆ > 0 does not require the calculation of the Hull-White extensions (θ (k) ) k∈N0 according to (10), but the knowledge of the parameter process (b(k), β(k), Σ(k)) k∈N0 is sufficient. The Hull-White extensions are fully encoded in the yield curve process (Y(k, ·)) k∈N0 , and we can avoid the inversion of (potentially) high dimensional matrices C(β(k)) k∈N0 .
Further remarks. • CRC of the multifactor Vasiček spot rate model can be defined directly in the HJM framework assuming a stochastic dynamics for the parameters. However, solely from the HJM representation, one cannot see that the yield curve dynamics is obtained, in our case, by combining well-understood Hull-White extended multifactor Vasiček spot rate models using the CRC algorithm of Section 3; that is, the Hull-White extended multifactor Vasiček model gives an explicit functional form to the HJM representation.
• The CRC algorithm of Section 3 does not rely directly on (ε * (t)) t∈N having independent and Gaussian components. The CRC algorithm is feasible as long as explicit formulas for ZCB prices in the Hull-White extended model are available. Therefore, one may replace the Gaussian innovations by other distributional assumptions, such as normal variance mixtures. This replacement is possible provided that conditional exponential moments can be calculated under the new innovation assumption. Under non-Gaussian innovations, it will no longer be the case that the HJM representation does not depend on the Hull-White extension θ (k) ∈ R N .
• Interpretation of the parameter processes will be given in Section 5, below.

Real World Dynamics and Market Price of Risk
All previous derivations were done under an equivalent martingale measure P * for the bank account numeraire. In order to statistically estimate parameters from market data, we need to specify a Girsanov transformation to the real-world measure, which is denoted by P. We present a specific change of measure, which provides tractable spot rate dynamics under P. Assume that (λ(k)) k∈N0 and (Λ(k)) k∈N0 are R n -and R n×n -valued F-adapted processes, respectively. Let (X (k)) k∈N0 be the factor process obtained by the CRC algorithm of Section 3.1. Then, we assume that the n-dimensional F-adapted process (λ(k) + Λ(k)X (k)) k∈N0 describes the market price of risk dynamics. We define the following P * -density process: (ξ(k)) k∈N0 The real-world probability measure P is then defined by the Radon-Nikodym derivative: An immediate consequence is that for k ∈ N 0 : has a standard Gaussian distribution under P, conditionally on F(k). This implies that under the realworld measure P, the factor process (X (k)) k∈N0 is described by: where we define: As in Assumption 1, we require Λ(k) to be such that the spectrum of α(k) is a subset of (−1, 1) n . Formula (14) describes the dynamics of the factor process (X (k)) k∈N0 obtained by the CRC algorithm of Section 3.1 under real-world measure P. The following corollary describes the yield curve dynamics obtained by the CRC algorithm under P, in analogy to Theorem 5.

Corollary 6.
Under real-world measure P satisfying (13), the yield curve dynamics (Y(k, ·)) k∈N0 obtained by the CRC algorithm of Section 3.1 has the following HJM representation for m > k + 1: , which are characterized by the market price of risk parameters λ(k) ∈ R n and Λ(k) ∈ R n×n .

Choice of Parameter Process
The yield curve dynamics obtained by the CRC algorithm of Section 3.1 require exogenous specification of the parameter process of the multifactor Vasiček Models (1) and (2) and the market price of risk process, i.e., we need to model the process: By Equation (9), the one-step ahead development of the CRC factor process X under P reads as: with F(t)-measurable parameters b(t), β(t) and Σ(t) and Hull-White extension θ (t) . Thus, on the one hand, the factor process (X (t)) t∈N0 evolves according to (17), and on the other hand, parameters (b(t), β(t), Σ(t), λ(t), Λ(t)) t∈N0 evolve according to the financial market conditions. Note that the process (θ (t) ) t∈N0 of Hull-White extensions is fully determined through CRC by (10). In order to distinguish the evolutions of (X (t)) t∈N0 and (b(t), β(t), Σ(t), λ(t), Λ(t)) t∈N0 , respectively, we assume that process (16) changes at a slower pace than the factor process, and therefore, parameters can be assumed to be constant over a short time window. This assumption motivates the following approach to specifying a model for process (16). For each time point t∆, we fit multifactor Vasiček Models (1) and (2) with fixed parameters (b, β, Σ, λ, Λ) on observations from a time window {t − K + 1, . . . , t} of length K. For estimation, we assume that we have yield curve observations ( y(k)) k=t−K+1,...,t = (( y 1 (k), . . . , y M (k))) k=t−K+1,...,t for times to maturity τ 1 ∆ < . . . < τ M ∆. Since yield curves are not necessarily observed on a regular time to the maturity grid, we introduce the indices τ 1 , . . . , τ M ∈ N to refer to the available times to maturity. Varying the time of estimation t∆, we obtain time series for the parameters from historical data. Finally, we fit a stochastic model to these time series. In the following, we discuss the interpretation of the parameters and present two different estimation procedures. The two procedures are combined to obtain a full specification of the model parameters.

Interpretation of Parameters
Thus, β determines the speed at which the factor process (X(t)) t∈N0 and the spot rate process (r(t)) t∈N0 return to their long-term means: A sensible choice of (β(t)) t∈N0 adapts the speed of mean reversion to the prevailing financial market conditions at each time point t∆.

Instantaneous
Variance. By Equation (3), we have under P * for t > 0: Thus, matrix Σ plays the role of the instantaneous covariance matrix of X, and it describes the instantaneous spot rate volatility.

State Space Modeling Approach
On each time window, we want to use yield curve observations to estimate the parameters of timehomogeneous Vasiček Models (1) and (2). In general, this model is not able to reproduce the yield curve observations exactly. One reason might be that the data are given in the form of parametrized yield curves, and the parametrization might not be compatible with the Vasiček model. For example, this is the case for the widely-used Svensson family [8]. Another reason might be that yield curve observations do not exactly represent risk-free zero-coupon bonds.
The discrepancy between the Vasiček model and the yield curve observations can be accounted for by adding a noise term to the Vasiček yield curves. This defines a state space model with the factor process as the hidden state variable. In this state space model, the parameters of the factor dynamics can be estimated using Kalman filter techniques in conjunction with maximum likelihood estimation ([9] Section 3.6.3). This is explained in detail in Sections 5.2.1-5.2.7 below.

Transition
System. The evolution of the unobservable process X under P is assumed to be given on time window {t − K + 1, . . . , t} by: with initial factor X(t − K) ∈ R n and parameters a = b − Σ 1 2 λ and α = β − Σ 1 2 Λ. The initial factor X(t − K) is updated according to the output of the Kalman filter for the previous time window {t − K, . . . , t − 1}. The initial factor is set to zero for the first time window available.
Remark. Parameters (b, β, Σ, λ, Λ) are assumed to be constant over the time window {t − K + 1,. . .,t}. Thus, we drop the index k compared to Equations (14) and (15). For estimation, we assume that the factor process evolves according to the time-homogeneous multifactor Vasiček Models (1) and (2) in that time window. The Hull-White extension is calibrated to the yield curve at time t∆ given the estimated parameter values of the time-homogeneous model.

Measurement System.
We assume that the observations in the state space model are given by: where: with A(·, ·) and B(·, ·) = (B 1 (·, ·), . . . , B n (·, ·)) given by Theorem 2 and M -dimensional F(k)-measurable noise term S 1 2 η(k) for non-singular S 1 2 ∈ R M ×M . We assume that η(k) is independent of F(k − 1) and ε(k) under P and that η(k) P ∼ N (0, 1). The error term S 1 2 η describes the discrepancy between the yield curve observations and the model. For S = 0, we would obtain a yield curve in (18) that corresponds exactly to the multifactor Vasiček one. Given the parameter and market price of risk value (b, β, Σ, λ, Λ), we estimate the factor using the following iterative procedure. For each fixed value of k ∈ {t − K, . . . , t} and fixed time t, we consider the σ-field F Y (k) = σ Y (s) t − K ≤ s ≤ k and describe the estimation procedure in this state space model.

Bayesian Inference in the Transition
System. The prediction error ζ(k) is used to update the unobservable factors.
where K(k) denotes the Kalman gain matrix given by: 5.2.6. Forecasting the Transition System. For the unobservable factor process, we have the following forecast: x 5.2.7. Likelihood Function. The Kalman filter procedure above allows one to infer factors X given the parameter and market price of risk values. Of course, in this section, we are interested in estimating these values in the first place. For this purpose, the procedure above can be used in conjunction with maximum likelihood estimation. For the underlying parameters Θ = (b, β, Σ, a, α), we have the following likelihood function given the observations ( y(k)) k=t−K+1,...,t : The maximum likelihood estimator (MLE) Θ is found by maximizing the likelihood function L t (Θ) over Θ, given the data. As in the EM (expectation maximization) algorithm, maximization of the likelihood function is alternated with Kalman filtering until convergence of the estimated parameters Θ MLE is achieved.

Estimation Motivated by Continuous Time Modeling
Rescaling the Time Grid. Assume factor process (X(t)) t∈N0 is given under P by X(0) ∈ R n and for t > 0: where a = b − Σ 1 2 λ and α = β − Σ 1 2 Λ. Furthermore, assume that α is a diagonalizable matrix with α = T DT −1 for T ∈ R n×n and diagonal matrix D ∈ (−1, 1) n×n . Then, the transformed process Z = (T −1 X(t)) t∈N0 evolves according to: where c = T −1 a and Ψ = T −1 Σ(T −1 ) . For d ∈ N + , the d-step ahead conditional distribution of Z under P is given by: s=0 D s ΨD s . Suppose we have estimated µ ∈ R n , the diagonal matrix γ ∈ (−1, 1) n and Γ ∈ R n×n on the time grid with size d∆, for instance, using MLE, as explained in Section 5.2. We are interested in recovering the parameters c, D and Ψ of the dynamics on the refined time grid with size ∆ from µ, γ and Γ. The diagonal matrix D and vector c are reconstructed from the diagonal matrix γ as follows: where logarithmic and power functions applied to diagonal matrices are defined on their diagonal elements. Note that for i, j = 1, . . . , n, we have: Therefore, we recover Ψ from γ and Γ as follows.
. From the formulas for c, D and Ψ, we observe that the F t−1 -conditional mean of D t Z: and the F t−1 -conditional volatility of D t Z: live on different scales as d → ∞; in fact, volatility dominates for large d. Under P for t > 0, we have: Therefore, setting D t X = X(t) − X(t − 1), we obtain as d → ∞:

Longitudinal Realized Covariations of Yields.
We consider the yield curve increments within the discrete-time multifactor Vasiček Models (1) and (2). The increments of the yield process (Y (t, t + τ )) t∈N0 for fixed time to maturity τ ∆ > 0 are given by: . For times to maturity τ 1 ∆, τ 2 ∆ > 0, we get under P: By Equation (20) for small grid size ∆, we estimate the last expression by: Formula (21) is interesting for the following reasons: • It does not depend on the unobservable factors X.
• It allows for direct cross-sectional estimation of β and Σ. That is, β and Σ can directly be estimated from market observations without knowing the market-price of risk.
• It is helpful to determine the number of factors needed to fit the model to market yield curve increments. This can be analyzed by principal component analysis.

•
It can also be interpreted as a small-noise approximation for noisy measurement systems of the form (18).
Let y 1 (k) and y 2 (k) be market observations for times to maturity τ 1 ∆ and τ 2 ∆ and at times k ∈ {t − K + 1, . . . , t}, also specified in Section 5.2. Then, the expectation on the left hand side of (21) can be estimated by the realized covariation: The quality of this estimator hinges on two crucial assumptions. First, higher order terms in (20) are negligible in comparison to Σ. Second, the noise term S 1 2 η in (18) leads to a negligible distortion in the sense that observations Y are reliable indicators for the underlying Vasiček yield curves.

Cross-Sectional
Estimation of β and Σ. Realized covariation estimator (22) can be used in conjunction with asymptotic relation (21) to estimate parameters β and Σ at time t∆ in the following way. For given symmetric weights w ij = w ji ≥ 0, we solve the least squares problem: where we optimize over β and Σ satisfying Assumption 1.

Inference on Market Price of Risk
Finally, we aim at determining parameters λ and Λ of the change of measure specified in Section 4. For this purpose, we combine MLE estimation (Section 5.2) with estimation from realized covariations of yields (Section 5.3). First, we estimate β and Σ by β RCov and Σ RCov as in Section 5.3. Second, we estimate a, b and α by maximizing the log-likelihood: for fixed β and Σ over b ∈ R n , a ∈ R n and α ∈ R n×n with spectrum in (−1, 1) The constraint on the matrix α ensures that the factor process is stationary under the real-world measure P. From Equation (15), we have λ = Σ − 1 2 (b − a) and Λ = Σ − 1 2 (β − α). This motivates the inference of λ by: and the inference of Λ by: We stress the importance of estimating as many parameters as possible from the realized covariations of yields prior to using maximum likelihood estimation. The MLE procedure of Section 5.2 is computationally intensive and generally does not work well to estimate volatility parameters.

Description and Selection of Data
We choose ∆ = 1/252, which corresponds to a daily time grid (assuming that a financial year has 252 business days). For the Swiss currency (CHF), we consider as yield observations the Swiss Average Rate (SAR), the London InterBank Offered Rate (LIBOR) and the Swiss Confederation Bond (SWCNB). See Figures 1 and 2.
• Short times to maturity. The SAR is an ongoing volume-weighted average rate calculated by the Swiss National Bank (SNB) based on repo transactions between financial institutions. It is used for short times to maturity of at most three months. For SAR, we have the Over-Night SARON that corresponds to a time to maturity of ∆ (one business day) and the SAR Tomorrow-Next (SARTN) for time to maturity 2∆ (two business days). The latter is not completely correct, because SARON is a collateral over-night rate and tomorrow-next is a call money rate for receiving money tomorrow, which has to be paid back the next business day. Moreover, we have the SAR for times to maturity of one week (SAR1W), two weeks (SAR2W), one month (SAR1M) and three months (SAR3M); see also [10].
• Short to medium times to maturity. The LIBOR reflects times to maturity, which correspond to one month (LIBOR1M), three months (LIBOR3M), six months (LIBOR6M) and 12 months (LIBOR12M) in the London interbank market.
These data are available from 8 December 1999, and we set 15 September 2014 to be the last observation date. Of course, SAR, LIBOR and SWCNB do not exactly model risk-free zero-coupon bonds, and these different classes of instruments are not completely consistent, because prices are determined slightly differently for each class. In particular, this can be seen during the 2008-2009 financial crisis. However, these data are in many cases the best approximation to CHF risk-free zero-coupon yields that is available. For the longest times to maturity of SWCNB, one may also raise issues about the liquidity of these instruments, because insurance companies typically run a buy-and-hold strategy for long-term bonds.
In Figures 3-6, we compute the realized volatility RCov(t, τ, τ ) 1 2 of yield curves ( y τ (k)) k=t−K+1,...,t for different times to maturity τ ∆ and window length K; see Equation (22). In Figures 2 and 6, we observe that SAR fits SWCNB better than LIBOR after the financial crisis of 2008. For this reason, we decide to drop LIBOR and build daily yield curves from SAR and SWCNB, only. The mismatch between LIBOR, SAR and SWCNB is attributable to differences in liquidity and the credit risk of the underlying instruments.
6.2.1. Discussion of Identification Assumptions. We select short times to maturity (SAR) to estimate parameters b, β, Σ, a and α. This is reasonable because these parameters describe the dynamics of the factor process and, thus, of the spot rate. As we are working on a small (daily) time grid, asymptotic Formulas (20) and (21) are expected to give good approximations. Additionally, it is reasonable to assume that the noise covariance matrix S in data-generating Model (18) is negligible compared to (21). Therefore, we can estimate the left hand side of (21) by the realized covariation of observed yields; see estimator (22). Then, we determine the Hull-White extension θ in order to match the prevailing yield curve interpolated from SAR and SWCNB.      Figure 6: A selection of SAR, LIBOR and SWCNB realized volatility RCov(t, τ, τ ) 1 2 for τ = 1, 63, 252, 504, window length K = 21 (lhs) and K = 126 (rhs). Note that LIBOR looks rather differently from SAR and SWCNB after the financial crisis of 2008.

Determination of the Number of Factors.
We need to determine the appropriate number of factors n. The more factors we use, the better we can fit the model to the data. However, the dimensionality of the estimation problem increases quadratically in the number of factors, and the model may become over-parametrized. Therefore, we look for a trade-off between the accuracy of the model and the number of parameters used. In Figure 7, we determine β 11 , . . . , β nn and Σ by solving optimization (23) numerically for three observation dates and n = 2, 3. A three-factor model is able to capture rather accurately the dependence on the time to maturity τ . In Figure 8, we compare the realized volatility of the numerical solution of (23) to the market realized volatility for all observation dates. We observe that in several periods, the two-factor model is not able to fit the SAR realized volatilities accurately for all times to maturities. The three-factor model achieves an accurate fit for most observation dates. The model exhibits small mismatches in 2001, 2008-2009 and 2011-2012. These are periods characterized by a sharp reduction in interest rates in response to financial crises. In September 2011, following strong appreciation of the Swiss Franc with respect to the Euro, the SNB pledged to no longer tolerate Euro-Franc exchange rates below the minimum rate of 1.20, effectively enforcing a currency floor for more than three years. As a consequence of the European sovereign debt crisis and the intervention of the SNB starting from 2011, we have a long period of very low (even negative) interest rates.      of X 2 and X 3 are higher than that of X 1 for most of the observation dates. We also see that the volatility of X 1 is lower than that of X 2 and X 3 . In 2011, we observe large spikes in the factor volatilities. Starting from 2011, we have a period with strong correlations among the factors. From these results, we conclude that the three-factor Vasiček model is reasonable for Swiss interest rates. Particularly challenging for the estimation is the period 2011-2014 of low interest rates following the European sovereign debt crisis and the SNB intervention. In Figure 9a (rhs), we observe that the difference in the speeds of meanreversion under the risk-neutral and real-world measures is negligible. The difference between b and a is considerable in certain time periods. From the estimation results, we conclude that a constant market price of risk assumption is reasonable and set from now on Λ = 0. In Figure 10, we compute the objective function of optimization (24) for (b, β, Σ, a, α) = (0, β RCov , Σ RCov , 0, β RCov ) and compare it to the numerical solution ( b MLE , β RCov , Σ RCov , a MLE , β RCov ). We observe that in 2003-2005 and 2010-2014, the parameter configuration (0, β RCov , Σ RCov , 0, β RCov ) is nearly optimal. In these periods, we have very low interest rates, and therefore, estimates of b and a close to zero are reasonable. Given the estimated parameters, we calibrate the Hull-White extension by equation (10) to the full yield curve interpolated from SAR and SWCNB; see Figure 11. We point out that our fitting method is not a purely statistical procedure; rather, it is a combination of estimation and calibration in accordance with the paradigm of robust calibration, as explained in [1].

Selection of a Model for the Vasiček Parameters.
In the following, we use the CRC approach to construct a modification of the Vasiček model with stochastic volatility. We model the process (Σ(t)) t∈N0 by a Heston-like [11] approach. We assume deterministic correlations among the factors and stochastic volatility given by:  where ϕ ∈ R 3 + , φ = diag (φ 11 , φ 22 , φ 33 ) ∈ R 3×3 , Φ 1 2 ∈ R 3×3 non-singular, and for each t ∈ N, ε(t) has a standard Gaussian distribution under P, conditionally given F(t − 1). Moreover, we assume that (ε(t), ε(t)) is multivariate Gaussian under P, conditionally given F(t − 1). Note that ε(t) and ε(t) are allowed to be correlated. The matrix valued process (Σ(t)) t∈N0 is constructed combining this stochastic volatility model with fixed correlation coefficients. This model is able to capture the stylized fact that volatility appears to be more noisy in high volatility periods; see Figure 9b. We use the volatility time series of Figure 9b to specify ϕ, φ and Φ. We rewrite the equation for the evolution of the volatility as: and use least square regression to estimate ϕ, φ and Φ. From the regression residuals, we estimate the correlations between ε(t) and ε(t). Figures 12-14 show the estimates of ϕ, φ and Φ.

Simulation and Back-Testing
Section 6.2 provides a full specification of the three-factor Vasiček CRC model under the risk-neutral and real-world probability measures. Various model quantities of interest in applications can then be calculated by simulation.
6.3.1. Simulation. The CRC approach has the remarkable property that yield curve increments can be simulated accurately and efficiently using Theorem 5 and Corollary 6. In contrast, spot rate models with stochastic volatility without CRC have serious computational drawbacks. In such models, the calculation of the prevailing yield curve for given state variables requires Monte Carlo simulation. Therefore, the simulation of future yield curves requires nested simulations.
6.3.2. Back-Testing. We backtest properties of the monthly returns of a buy and hold portfolio investing equal proportions of wealth in the zero-coupon bonds with times to maturity of 2, 3, 4, 5, 6 and 9   Figure 14: Estimation of correlations ρ 21 , ρ 31 and ρ 32 (lhs) and correlations Cor [ε(t), ε(t) | F(t − 1)] (rhs). We use a time window of 252 observation for the regression. The residuals ε are calculated using the parameter estimates of Figure 9. months and 1, 2, 3, 5, 7 and 10 years. We divide the sample into disjoint monthly periods and calculate the monthly return of this portfolio assuming that at the beginning of each period, we invest in the bonds with these times to maturity in equal proportions of wealth. The returns and some summary statistics are shown in Figure 15. We observe that the returns are positively skewed, leptokurtic and have heavier tails than the Gaussian distribution. These stylized facts are essential in applications.
For each monthly period, we select a three-factor Vasiček model and its CRC counterpart with stochastic volatility. Then, we simulate for each period realizations of the returns of the test portfolio. By construction, the Vasiček model generates Gaussian log-returns and is unable to reproduce the stylized facts of the sample; see Tables 1 and 2 and Figure 16. Increasing the number of factors does not help much, because the log-returns remain Gaussian. On the other hand, CRC of the Vasiček model with stochastic volatility provides additional modeling flexibility. In particular, we can see from the statistics in Table 2 and the confidence intervals in Figure 16 that the model matches the return distribution better than the Vasiček model. As explained in Figure 16, statistical tests assuming the independence of disjoint monthly periods show that the difference between the Vasiček model and its CRC counterpart is statistically significant. We conclude that the three-factor CRC Vasiček model is a parsimonious and tractable alternative that provides reasonable results.
Vasicek   Monthly portfolio logarithmic returns Figure 15: Logarithmic monthly returns of a buy and hold portfolio investing in equal wealth proportions in the zero-coupon bonds with times to maturity of 2, 3, 4, 5, 6 and 9 months and 1, 2, 3, 5, 7 and 10 years. For each monthly period, we calculate the logarithmic return of this portfolio assuming that at the beginning of each period, we are invested in the bonds with these times to maturity in equal proportions of wealth.

Regulatory Framework.
The type of analysis that was performed in the previous section is an integral component of the present regulatory framework for risk management. In the Basel framework [5], the capital charge for the trading book is based on quantile risk measures. Under the internal model approach ( [5], Section 2.VI.D), a bank calculates quantiles for the distribution of possible 10day losses based on recent market data under the assumption that the trading book portfolio is held fixed over the time period. The approach relies on accurate modeling of the distribution of portfolio returns over holding periods of multiple days. A similar analysis is required by the Basel ( [5], Section 2.VI.D) regulatory framework for model validation and stress testing: model validation is performed by backtesting the historical performance of the model, and stress tests are carried out using the same methodology by calibrating the model to historical periods of significant financial stress. These tasks can be accomplished using the CRC approach by selecting suitable classes of affine models and parameter processes. The approach is fairly general, since there are few restrictions on the parameter processes. In particular, it allows for stochastic volatility and can be used to create realistic non-Gaussian distributions of multi-period bond returns (see Section 6.3.2). Nevertheless, computing these bond return distributions does not require nested simulations. This is crucial for reasons of efficiency. Moreover, the flexibility in the specification of the parameter processes makes the CRC approach well suited for stress testing, because it allows one to freely select and specify stress scenarios.

Conclusions
• Flexibility and tractability. Consistent re-calibration of the multifactor Vasiček model provides a tractable extension that allows parameters to follow stochastic processes. The additional flexibility can lead to better fits of yield curve dynamics and return distributions, as we demonstrated in our numerical example. Nevertheless, the model remains tractable. In particular, yield curves can be simulated efficiently using Theorem 5 and Corollary 6. This allows one to efficiently calculate model quantities of interest in risk management, forecasting and pricing.
• Model selection. CRC models are selected from the data in accordance with the robust calibration principle of [1]. First, historical parameters, market prices of risk and Hull-White extensions are inferred using a combination of volatility estimation, MLE and calibration to the prevailing yield curve via Formulas (23-26, 10). The only choices in this inference procedure are the number of factors of the Vasiček model and the window length K. Then, as a second step, the time series of estimated historical parameters are used to select a model for the parameter evolution. This results in a complete specification of the CRC model under the real world and the pricing measure.
• Application to modeling of Swiss interest rates. We fitted a three-factor Vasiček CRC model with stochastic volatility to Swiss interest rate data. The model achieves a reasonably good fit in most time periods. The tractability of CRC allowed us to compute several model quantities by simulation. We looked at the historical performance of a representative buy and hold portfolio of Swiss bonds and concluded that a multifactor Vasiček model is unable to describe the returns of this portfolio accurately. In contrast, the CRC version of the model provides the necessary flexibility for a good fit.
Proof of Theorem 5. We add and subtract −A (k) (k, m) + B (k) (k, m) X (k) to the right hand side of Equation (12) and obtain: We have the following two identities from Section 3.1.