{'fraction_non_alphanumeric': 0.029723403719829865, 'fraction_numerical': 0.011212342385422696, 'mean_word_length': 4.4578983516483515, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Error Analysis of a Model Order Reduction Framework for Financial Risk Analysis

A parametric model order reduction (MOR) approach for simulating the high dimensional models arising in financial risk analysis is proposed on the basis of the proper orthogonal decomposition (POD) approach to generate small model approximations for the high dimensional parametric convection-diffusion reaction partial differential equations (PDE). The proposed technique uses an adaptive greedy sampling approach based on surrogate modeling to efficiently locate the most relevant training parameters, thus generating the optimal reduced basis. The best suitable reduced model is procured such that the total error is less than a user-defined tolerance. The three major errors considered are the discretization error associated with the full model obtained by discretizing the PDE, the model order reduction error, and the parameter sampling error. The developed technique is analyzed, implemented, and tested on industrial data of a puttable steepener under the two-factor Hull-White model. The results illustrate that the reduced model provides a significant speedup with excellent accuracy over a full model approach, demonstrating its potential applications in the historical or Monte Carlo value at risk calculations.


Introduction
Investors use risk analysis to determine whether to undertake a particular venture, the plausible return and how to mitigate potential losses of an activity. The risk analysis of financial instruments often requires the valuation of such instruments under a wide range of future market scenarios. A suitable valuation function takes these market scenarios as input parameters and determines the fair value of financial instruments. Examples of such risk analysis tasks are the calculation of Value-at-Risk (VaR) or the expected shortfall to estimate worst-case scenarios of financial holdings. Package retail investments and insurance products (PRIIPs) are the packaged instruments that are offered to retail investors. In order to make PRIIPs from different manufacturers more comparable concerning their risks and returns, the European regulation (EU) 1286/2014 requires manufacturers of PRIIPs to supply key information documents (KIDs) to possible retail investors that are easy to read and to understand [17]. The commission delegated regulation (EU) 2017/653 formulates the details of how the risk and the possible returns of a PRIIP have to be calculated [18]. There are four different categories of PRIIPs. In this paper, we concentrate on category 3 instruments, for which at least 10 000 market data scenarios have to be generated, and the PRIIP has to be valuated under these scenarios.
In our approach, the financial instruments are evaluated via the dynamics of short-rate models [11], based on convection-diffusion-reaction partial differential equations (PDE). The choice of the short-rate model depends on the underlying financial instrument. Some of the prominent financial models are the one-factor Hull-White model [26], the shifted Black-Karasinski model [8], and the two-factor Hull-White model [28]. These models are calibrated based on market data like yield curves that generate a high dimensional parameter space [11]. To perform the risk analysis, the financial model needs to be solved for such a high dimensional parameter space, and this requires efficient algorithms. For instruments with high complexity and long-time horizons, computing times of minutes for a single valuation are not unusual.
Thus, in a previous paper [7], we have established a parametric model order reduction approach based on a variant of the proper orthogonal decomposition approach, which significantly reduces the overall computation time [6,13]. This nonlinear MOR approach is computationally feasible [23] as it determines low dimensional linear (or affine) subspaces [44,57] via a truncated singular value decomposition of a snapshot matrix [53], that is computed by simulating the full model obtained by discretizing the PDE for a small number of pre-selected training parameter values. The question of how to select these parameters is often the most challenging part of the process. Our previous work [7] has established greedy algorithms to determine the training parameters more efficiently. The adaptive greedy approach searches the entire parameter space efficiently using a surrogate model and determines the best suitable training parameters. The principal component regression technique has been used to construct the surrogate model. This paper provides a detailed error analysis of the model order reduction framework. We have considered the hierarchical framework to perform the error analysis that starts by considering the solution of the partial differential equation as the true solution. We then discretize the underlying PDE to generate the full model and project it on a reduced basis to obtain the reduced model. Each stage described here is associated with some numerical error. Three major numerical errors considered are the discretization error associated with the full model, the model order reduction error, and the sampling error associated with the sampling algorithms. The total error associated with the model order reduction framework is the sum of these three significant errors. Our approach aims to minimize the total error such that it is less than the user-defined tolerance generating the suitable reduced model. We define a discretization error estimator based on estimating the exact solution of the PDE, which is higher-order accurate than the underlying numerical solution [48]. The method is known as Richardson extrapolation [45], which uses a sequence of solutions to estimate the discretization error and can be applied to almost any type of discretization approach as well as to both local and global quantities [49]. The discretization error is mainly governed by the choice of the grid size, so we present an approach to select an optimal grid based on the developed error estimator.
The model order reduction error is composed of the projection error associated with the proper orthogonal decomposition approach and how well the reduced model approximates the full model [33]. Both the projection error and the reduced model error depend on the dimension of the reduced basis. We present an algorithm to select an optimal reduced dimension that minimizes the projection error and the reduced model error. To obtain a reduced basis, we have to compute a singular value decomposition of the snapshot matrix, which is a computationally costly task. Recent research has shown that the randomized singular value decomposition is a good alternative and it is noticeably faster than the basic SVD for some applications [36,42,56]. Thus, we use the randomized truncated singular value decomposition, developed in [24], in the adaptive greedy sampling approach to determine the reduced basis. The sampling algorithm used to obtain the training parameters to construct a reduced basis introduces a sampling error. In the adaptive greedy sampling, the projection error due to the principal component analysis (PCA) and the prediction error of the surrogate model form the total sampling error. We analyze the projection error due to the PCA in a similar sense to the projection error of the POD, while the mean square error of prediction is used to assess the quality of the surrogate model.
We also perform the sensitivity analysis for these errors to know the relative contribution of each error to the total error [51]. It provides useful insight into which error contributes most to the variability of the model output [14]. The sensitivity analysis assists in improving the model output by selecting the proper grid size, reduced dimension, and training parameters. It also helps to allocate computational resources efficiently as more computational effort could be spent on the most sensitive parameters.
In a previous paper, we have considered a numerical example of a floater with caps and floors under the one-factor Hull-White model [7]. However, there are more complicated instruments, e.g., steepeners whose coupons depend on the difference between the two interest rates [1]. For such instruments, we need to consider a multi-factor model. We solve a numerical example of the puttable steepener with caps and floors using the two-factor Hull-White model in this paper.

Model Order Reduction
In an affine one-factor model, the forward rates f (t, T 1 ), f (t, T 2 ) at two different tenors T 1 , T 2 are perfectly correlated, i.e., their correlation coefficient is one, which is unrealistic [9]. This perfect correlation only allows parallel shifts of the yield curves. However, we often noticed the yield curve steeping where short-term rates are low and the long-term rates are high [11]. Also, for instruments like steepeners [2], where coupons depend on the difference between two rates, we need to consider a multi-factor model like the two-factor Hull-White model. The two-factor model overcomes these drawbacks by adding a stochastic disturbance to the drift term. In the two-factor Hull-White model, the short rate r is assumed to satisfy the following stochastic differential equation [27] dr(t) = (θ(t) + u(t) − α(t)r(t))dt + σ 1 (t)dW 1 (t), where α, b, σ 1 , σ 2 > 0 and W 1 , W 2 are two Brownian motions under the risk neutral measure such that where −1 ≤ γ ≤ 1 is the correlation coefficient. To evaluate the underlying financial instrument, we use the partial differential equation for the two-factor Hull-White model and then solve it using appropriate boundary conditions and known terminal conditions. Consider a financial instrument V (t, r(t), u(t)) contingent on the stochastic interest rate movement r(t).
The two-factor Hull-White PDE is then given as [27] ∂V ∂t (2) and the deterministic function θ(t) is chosen to fit the simulated yield curves. However, the regulation does not provide a methodology for obtaining the parameters α(t), b(t), σ 1 (t), σ 2 (t). Thus, considering them as constants should lead to more robust results than time-dependent parameters [1]. We consider the resulting robust two-factor Hull-White model with parameters α, b, σ 1 , σ 2 > 0 as positive constants and only the parameter θ(t) as time-dependent. Our results can, however, be extended to the more general case. According to the PRIIPs regulations, we have to perform at least s = 10 000 yield curve simulations. We construct a simulated yield curve matrix which is then used to calibrate the parameter θ(t). The calibration based on s = 10 000 different simulated yield curves generates s different piecewise constant parameters θ ℓ (t), which change their values θ ℓ,i only at the m tenor points. We incorporate these in a matrix A detailed procedure for the yield curve simulation and the calibration of θ(t) based on simulated yield curves is described in [7,38], respectively. We have applied a finite element method to solve the PDE numerically, for more details see [1]. This discretization is a parametric high dimensional model of the form with given terminal vector V T , and matrices A(ρ ℓ ) ∈ R M ×M , and B(ρ ℓ ) ∈ R M ×M . Here ρ = {α, b, σ 1 , σ 2 , γ, θ(t)} is the group of model parameters. We call this the full model (FM) for the model reduction procedure. We solve (3) by propagating backward in time. Here again ℓ = 1, . . . , s = 10 000, m is the total number of tenor points, and we need to solve this system at each time step n with an appropriate boundary conditions and a known terminal value for the underlying instrument. Altogether we have a parameter space P of size 10 000 × m to which we now apply model reduction.
To perform the parametric model reduction for system (3), we employ Galerkin projection onto a low dimensional subspace viaV where the columns of Q ∈ R M ×d represent the reduced basis with d ≪ M , V n d is a vector of reduced coordinates, andV n ∈ R M is the solution in the nth time step obtained using the reduced model. For the Galerkin projection we require that the residual of the reduced state is orthogonal to the reduced basis matrix Q, i.e., so that by multiplying p n (V n d , ρ ℓ ) with Q T , we get where A d (ρ ℓ ) ∈ R d×d and B d (ρ ℓ ) ∈ R d×d are the parameter dependent reduced matrices. We obtain the Galerkin projection matrix Q in (4) based on a proper orthogonal decomposition (POD) approach, which generates an optimal order orthonormal basis Q in the least square sense that is independent of the parameter space P and we do this by the method of snapshots. The snapshots are nothing but the state solutions obtained by simulating the full model for selected parameter groups. We assume that we have a training set of parameter groups ρ 1 , . . . , ρ k ∈ [ρ 1 , ρ s ]. We compute the solutions of the full model for this training set and combine them in a snapshot matrixV = [V (ρ 1 ), V (ρ 2 ), ..., V (ρ k )]. The POD method solves, see [22], for an orthogonal matrix Q ∈ R M ×d via a truncated singular value decomposition (SVD) where φ i and ψ i are the left and right singular vectors of the matrixV respectively, and Σ i are the singular values. The truncated SVD computes only the first k columns of the matrix Φ. The detailed error analysis of the model order reduction and the randomized truncated singular value decomposition is explained in Subsection 3.2. We choose only d out of k POD modes to construct Q = [φ 1 · · · φ d ] which minimizes the projection error (28) We summarize the procedure of selecting the dimension d of the reduced basis in Algorithm 5. It is evident that the quality of the reduced model strongly depends on the selection of parameter groups ρ 1 , ..., ρ k that are used to compute the snapshots. Hence it is essential to introduce an efficient sampling technique for the parameter space. We could consider standard sampling techniques, like uniform sampling or random sampling [31]. However, these techniques may neglect vital regions within the parameter space. As an alternative, a greedy sampling method has been suggested in the framework of model order reduction. The greedy sampling technique selects the parameters that maximize the error between the reduced model and the full model. Further, the reduced basis is obtained using these selected training parameters. The relative error calculation is expensive as it requires the full model solution, so instead, the residual error associated with the reduced model (5) is used as an error estimator. However, it is not reasonable to compute an error estimator for the entire parameter space. This means that one has to select a pre-defined parameter set as a subset of the high dimensional parameter space to train the greedy sampling algorithm. We usually select this pre-defined subset randomly. A random selection may neglect the crucial parameters within the parameter space. Thus, to surmount this problem, we implemented an adaptive greedy sampling approach. The algorithm chooses the most suitable parameters adaptively using an optimized search based on surrogate modeling. This approach evades the cost of computing the error estimator for each parameter within the parameter space and instead uses a surrogate model to locate the training parameter set. We have built the surrogate model using the principal component regression technique. See [7] for detailed sampling algorithms and their implementation. This paper aims to study the numerical errors that arise during different stages of the model order reduction approach.

Error Analysis
In most cases, the exact solutions of mathematical models are unknown, and they incorporate some modeling errors. We have no other way to compute approximate solutions and to analyze them. However, the approximate solution contains some numerical errors. It is necessary to reduce these numerical errors as much as possible to improve the level of accuracy of the numerical solutions. Our mathematical model is a partial differential equation (PDE) with appropriate boundary conditions and a known terminal condition. To solve this system, we have designed an approximate solution by discretizing the PDE and thus committing a discretization error. Using model reduction based on the POD leads to a projection error, which depends on the excluded singular values while generating the reduced basis. To know how well the reduced model approximates the full model, we have to determine the model order reduction error. In this work, the training parameters are chosen based on either the classical greedy sampling or the adaptive greedy sampling approaches. Thus, we also have to consider a sampling error associated with the sampling algorithms as well. Figure 1 shows the model hierarchy to obtain the reduced model with associated errors in each stage. Usually, the formulated mathematical model is not an exact representation of the underlying phenomenon but subject to many simplifying assumptions. In this paper, we assume the solution of the partial differential equation as the true solution, i.e., the modeling error ǫ m will be neglected. Let Furthermore, the reduced model M d h of dimension d is obtained. Let theV be the result obtained from the reduced model. The reduced model error can then be given as The result obtained using the reduced model encloses the total numerical error, and based on that, we can write the final output equation with numerical errors The reduced model is obtained based on the proper orthogonal decomposition approach, which induces a projection error ǫ POD . The total model order reduction error ǫ d is then the sum of a projection error ǫ POD , and a reduced model error ǫ RM We use a sampling approach to select the most relevant parameters. This produces a reasonably accurate reduced basis, but the sampling algorithms contribute to an additional sampling error. Thus, the total simulation error ǫ T can be estimated by the discretization error ǫ h , a sampling error ǫ samp and the model order reduction error ǫ d as We aim to generate the reduced model and to solve the financial instrument for the entire parameter space such that the as a cost function the total error is minimized. This total error depends on the grid size h, the reduced dimension d, and the quality of training parameters c used to obtain the reduced basis. We optimize these variables such that the total error ǫ T is less than the user-defined tolerance e tol , i.e., min {h,d,c} Another important issue is to know the contribution and the effect of each defined error on the total error (11). Sensitivity analysis enables identifying the errors that have the most significant influence on the model output. The sensitivity approach presented in Subsection 3.4 ranks the numerical errors according to their contribution to the total error and helps to efficiently allocate the computational resources. The extra computational effort can be spent on the cost parameters most sensitive to the final output at the expense of the computational cost required to obtain the least sensitive parameters.

Discretization of a PDE
For the discretization of a PDE, we use a finite dimension grid sd k=1 [r k , u k ] × [0, T ] in sd spatial dimensions and one temporal dimension. The simulation time [0, T ] is divided into N equal intervals, [t n , t n−1 ] where n = 1, · · · , N . We have considered an implicit time-stepping scheme. To obtain an optimal solution of the discrete model, one must determine the best suitable time step and the grid size. We implement an adaptive time-stepping scheme to achieve the optimal time step, which solves the model at the desired accuracy but with low computational cost while the grid size is selected such that the discretization error is less than a specified user-defined tolerance. We have a PDE of the form where L(r, u) is a linear differential operator containing all spatial derivatives. The idea is to obtain the optimal time step ∆t opt such that it minimizes the relative error between the exact solution V and the solution obtained using the time step ∆t opt Assuming the numerical method used is of order p in time, the exact value of the PDE V , and the approximated value V h (∆t) with the time step ∆t, the error for the approximation is given as [55] By neglecting the higher order terms we can approximately determine the value of the constant C as Let V h (K∆t) be the solution obtained using the time step K∆t Considering a time discretization scheme of order p = 1 and substituting (14) into (13), we get We calculate ∆t opt such that the relative error ǫ t is less than the user-defined tolerance where . Thus, we get the relative error estimateǫ t by substituting C as To determine the optimal step size we setǫ which yields Since the roundoff or iterative convergence errors that arise while solving the system of linear equations are negligible compared to the discretization error, we can neglect them [35]. The adaptive time-stepping strategy aims to find the time step that will achieve the desired accuracy at a relatively low cost [5,16]. Thus, we have to determine the largest possible time step that satisfies (16). An automatic time step control can be executed using the optimal time step ∆t opt , as shown in Algorithm 1. The step size controller is designed on the assumption that the largest possible step size should be selected. If ∆t opt > ∆t, we can say that we have a good estimation of the solution that satisfies the user-defined tolerance. Thus, the algorithm runs until the optimal step size ∆t opt becomes larger than the step size ∆t.  Compute an optimal time step ∆t opt = e t tol ∆t 2 (K 2 −1) Increment K = round(∆t/∆t opt ) + 1.
It is important to note that we have to select the time step such that all key time points are achievable, i.e., coupon dates, put dates, or valuation dates.
We also need a discretization error estimator based on estimating the exact solution for the differential equation, see [48]. We use the Richardson extrapolation approach that estimates the exact solution using the formal rate of convergence of a discretization method and two or three different solutions on systematically refined meshes [49]. One would argue that the computation of multiple simulations in the asymptotic range of the full model increases the computational burden. However, the main advantage of the Richardson approach is that it can be used as a post-processing technique and can be applied to any discretization method (e.g., finite element or finite difference). It estimates the total error, including locally generated errors and those transported from other regions of the domain, see also [3,34,43].
Consider a discretization error as where V is the exact solution of the mathematical model and V h is a solution obtained on a grid spacing h using a discretized model. In the case of a two-dimensional system where E is the number of elements, and A (e) i is the area of ith element. For a smooth enough solution with no discontinuities, we can expand the numerical solution V h using Taylor expansion which leads to the discretization error by neglecting higher order terms. Here p f is the formal order of accuracy determined by the preferred discretization scheme and is the first term in the series. To accurately estimate the discretization errorǭ h , Richardson extrapolation demands that all other sources of numerical errors (e.g., round-off error, iterative convergence errors) are negligible compared to the discretization error [35]. Consider a grid refinement index g as the ratio of coarse to fine mesh size. It typically is recommended see e.g. [12], to use the grid refinement ratio between 1.5 to 2 as it allows to distinguish the discretization error from other error sources.
Let h fine = h, and h coarse = gh, then the discretization error estimator is [48] To access the confidence in the discretization error estimate, it is necessary to define the observed order of accuracyp [48]. The requirement of the reliability of the discretization error estimate is that the solutions are in the range of asymptotic convergence. To demonstrate that the asymptotic range for different solutions has been achieved, we must have at least three different discretized solutions with three different mesh sizes. We further calculate the observed order using these discretized solutions over a range of meshes. According to [12], when the observed order of accuracy is equal to the formal order of accuracy, we can have a high degree of confidence in the error estimator. We can determine the observed order of accuracy as follows. Consider three different mesh sizes h 1 , h 2 and h 3 such that Let h 1 = h then h 2 = g 12 h and h 3 = g 23 g 12 h. Using the formal power series and subtracting (20) from (21) and (19) from (20) and settingp ≈ p, we obtain Dividing (22) by (23) and neglecting higher order terms we have We can then determinep by an iterative procesŝ where an initial guess forp is the formal order of accuracy. Algorithm 2 shows the procedure to compute the observed order of accuracy. We consider it converged when the change in two subsequentp k+1 andp k is not significant.
Algorithm 2 (Algorithm to estimate the observed order of accuracy).
Substituting p =p and g = g 12 in (18) leads to a discretization error estimate [48] ǫ h = 1 When the observed order of accuracyp matches the formal order of accuracy p f , we have high confidence that this error estimate is accurate. However, it may happen that the observed order of accuracy does not match the formal order of accuracy [43], then the error estimator is not reliable. This may be due to the failure of obtaining the solutions within the asymptotic range or the lack of additional information. Also, the implemented boundary conditions and the terminal/initial conditions change the observed order of accuracy. However, we can fix this problem by converting the error estimator into an epistemic uncertainty estimator [48]. Epistemic uncertainties are different from random or aleatory uncertainties, as epistemic uncertainties arise due to the lack of knowledge. A grid convergence index (GCI) measures how far the computed value is away from the asymptotic numerical value. The GCI converts the error estimate into an error or uncertainty band, which is appropriate when one does not have a high degree of confidence in the error estimate. A small value of the GCI indicates that the computation is within the asymptotic range [47]. The GCI for the fine grid solution is given as [46] where • F s = 3 if the observed order of accuracy is calculated using two different solutions.
• F s = 1.25 if the observed order of accuracy is calculated using three different solutions and when the formal order agrees with the observed order within 10%.
We cannot ensure that the approximate solution via Richardson extrapolation approximates the exact solution of mathematical model well. Thus, it is beneficial to use a safety factor to avoid the failure of the designed error estimator. For the first case, one cannot guaranty that the solution obtained with only two grid sizes is within the asymptotic range. Thus, it is recommended to use the GCI with caution. In [48] a safety factor of 3 is suggested. In the second case, if the formal order of accuracy matches the observed order of accuracy, then it is recommended to use a safety factor of 1.5. In this paper, whenp agrees p f within 10%, then the factor of safety of 1.25 is used. Algorithm 3 shows the steps to obtain the suitable full model dimension. The algorithm is initiated by reasonably moderate grid size h init 1 . We construct two new grid sizes h 2 , h 3 based on the grid refinement ratios g 12 , g 23 . Furthermore, a full model is solved for h 1 , h 2 , h 3 to determine the observed order of accuracyp. Once we calculatep, the algorithm computes the discretization error estimatorǭ h and checks it against the user-defined tolerance e h tol . The algorithm is terminated ifǭ h < e h tol else, it refines h 1 , and the process repeats itself until convergence or after K iterations, and the last grid size is used to construct the full model.
Solve the full model with a grid size h 1 and store the result in V 1 . Compute Solve the full model with grid sizes h 2 and h 3 and store the results in V 2 and V 3 respectively.
Compute the discretization error estimator given by (25)

END.
Refine the mesh and choose a finer grid size h 1 .

Model Order Reduction Error
For the model redeuction procedure we use a snapshot-based POD method, which represents the PDE solution in the space spanned by a reduced basis. The reduced basis is obtained by a linear combination of the snapshots generated by solving the full model for some training parameters.
Let X be the Hibert space given with an inner product (·, ·) X and a norm · X . Let V 1 , V 2 , · · · , V k ∈ X be the snapshots obtained by solving the full model. We consider as a span of the snapshots {V j } k j=1 , provided at least one of them is non-zero. Let denote an orthonormal basis of V with d = dim(V). The goal of POD is to identify the most important characteristics of an ensemble of snapshots. We can represent each (block) column of the snapshot matrix as The POD approach obtains a reduced basis such that for every ℓ ∈ {1, · · · , d}, the mean square error between the snapshot V j and the corresponding lth partial sum of (27) is minimized on average via is obtained by computing a truncated singular value decomposition of the snapshot matrixV [22] Let Σ 1 ≥ Σ 2 ≥ · · · Σ k > 0 denote the singular values, and φ 1 , φ 2 , · · · , φ k ∈ R n the left singular vectors ofV . The reduced basis Q then consists of d left singular vectors where the dimension d is chosen such that it minimizes the projection error ǫ POD [33] This error characterizes the ability of the snapshot data to be represented in a low dimensional space. The goal of reduced modeling is to select the best suitable snapshots and to reduce the POD error. We use the error estimator ǫ POD to choose the dimension d of the reduced basis in Algorithm 5. Furthermore, in order to test how well the reduced model approximates the full model, we define an error between the full model and the reduced model where V is the solution obtained using the full model andV is the result obtained using the reduced model. The quality of the reduced basis Q depends onV = [V 1 (ρ 1 ), v 2 (ρ 2 ), · · · , V k (ρ k )], which depend on the training parameter set {ρ 1 , ρ 2 , · · · , ρ k }. In [7], a classical and an adaptive greedy sampling approach is presented, to locate the training parameter set efficiently. We use an a posteriori residual error estimator (5) ε = p(·, ρ) to monitor the convergence of the classical greedy sampling algorithm and to terminate the training procedure. In the adaptive greedy sampling approach, the algorithm constructs an error modelǭ RM based on the available data E p from the residual errors {ε bef ,aft i } Imax i=1 and the reduced model errors {ǫ bef,aft RM,i } Imax i=1 and uses this error model to monitor the convergence of the algorithm. Here the superscripts bef and aft show the error values obtained before and after updating the reduced basis.
For this approach, one full model is solved for an optimal parameter group ρ I at each iteration i of the classical and the adaptive greedy procedures. The solution obtained for ρ I is then stored as V i (ρ I ) and updates the snapshot where T is the maturity. Now we consider as a span of these snapshots and obtain the reduced basis by computing its truncated singular value decomposition. At each greedy iteration, a new reduced basis Q i is obtained of dimension d and the projection error based on (28) is which is the sum of squares of neglected singular values. The computation of the SVD of an m × n matrix has an overall computational complexity of O(mn 2 ) [21], so for the matrixV ∈ R M ×(T ×i) it is After i iterations, the computational complexity due to the series of SVDs for the snapshot matrix will be We see that the complexity increases substantially with each proceeding iteration, which motivates us to implement a truncated SVD approach. Recent studies have shown that the randomized algorithms in some applications [19,24] yield an incredible speed-up. It is proposed to compute a truncated singular value decomposition using randomized algorithms in the following two steps. The approach starts by computing an approximate basis G ∈ R M ×k with k orthonormal columns for the range ofV such that where k < T × i and ǫ tol is some specified tolerance. In the second step the matrix G is used to compute the truncated SVD of the matrixV . LetṼ = G * V ∈ R k×(T ×i) be the reduced matrix determined fromV . Now we compute the SVD of the reduced matrix which leads to a substantial speedup. Multiplying both sides by G, we get Considering Φ = GΦ, we get our left singular vector matrix, which we can use to obtain a reduced basis. We can achieve even more speedup by incorporating the economy-sized (compact) SVD of the matrixṼ , which removes extra rows or columns of zeros from the diagonal matrix of singular values and the columns in eitherΦ orΨ. Removing these zeros and columns can improve execution time and reduces storage requirements without compromising the accuracy of the decomposition. The SVD determines the rank k of a matrix, as the number of nonzero singular values of a matrixṼ ∈ R k×T ×i . Then its compact SVD reads The left and right singular vector matrices areΦ ∈ R k×k , andΨ ∈ R (T ×i)×k . Thus, the desired left singular matrix is Φ = GΦ ∈ R M ×k , which we use to obtain the reduced basis Q. We obtain the orthonormal matrix G with as few columns as possible such that where ε tol is some tolerance and · is the L 2 -norm. The range of G is a k dimensional subspace that represents most of the action ofV , and we select k as small as possible. We can use the SVD to decide the dimension k. Let Σ j be the jth singular value of the matrix V . For each j ≥ 0 min where the columns of G are k dominant singular vectors ofV . The number of singular values ofV that exceeds the tolerance ε tol makes the minimal rank k. Draw k Gaussian random vectorsĜ = [ĝ 1 , · · · ,ĝ k ] ∈ R (T ×i)×k and project them on to the linear map of a matrixV such that We find the orthonormal basis of the matrix Y ∈ R M ×k , which is nothing but our desired matrix G, using a Gram-Schmidt like algorithm or QR factorization. We can determine G adaptively using the error (I − GG * )V . We can obtain some information about this error by calculating (I − GG * )Vĝ . For k Gaussian random vectors, [24] gives an posterior error estimator as This error estimate can be combined with any method for constructing an approximate basis for the range of Y , see [24] for more details. The adaptive algorithm 4 combined with Gram-Schmidt orthonormalization is given below. The algorithm is initiated by drawing k Gaussian random vectors of length T × i and computes the orthonormal set of vectors {g 1 , ..., g k } that form the matrix G. We stop the procedure once the largest value in a set { y i } i+k i is less than ǫ/(10 2/π).
IF max{ y i , y i+1 . . . , y i+k } < ǫ/(10 2/π) We use the randomized singular value decomposition in the greedy sampling approach to determine the reduced basis, as shown in Algorithm 5. The algorithm constructs a snapshot matrixV composed of full model solutions for each optimal parameter group at the ith iteration. It then selects k Gaussian random vectors and obtains the matrix G, using the function GSorth().
Algorithm 5 (Algorithm to obtain a reduced basis using a randomized singular value decomposition). Input: Snapshot matrixV , e d tol . Output: Q.
Obtain G using 'GSorth' function defined in Algorithm 4.
Choose a parameter group ρ test for testing.
Solve the full model for ρ test and store result in V .
Solve the reduced model for ρ test for store result inV .
A low dimensional matrixṼ is constructed using the newly obtained G with only k columns. The algorithm computes a compact SVD ofṼ and sets Φ = GΦ to obtain a reduced basis Q. Finally, the dimension of a reduced basis is selected that minimizes the projection error and the reduced model error.

Sampling Error
To design a surrogate model to locate the training parameters we will use principal component regression. There are two errors associated with the principal component regression model, i.e., a projection error due to the principal component analysis (PCA) and an error of prediction for the surrogate model.
Let us first discuss the projection error. PCA reduces the problem of estimating m coefficients to the more straightforward problem of determining p coefficients by selecting a few relevant principal components. It is a technique for projecting multidimensional data onto a low dimensional subspace with minimal loss of variance. We determine the projection error associated with the PCA in a similar sense to the POD. At the kth iteration we have a data matrixP k = [ρ 1 , . . . , ρ c k ] composed of c k parameter groups. To build the surrogate model, we compute an SVD of the matrixP k P k =ΦΣΨ, where the columns ofP kΨ are principal components. We use only p principal components to construct a fairly accurate surrogate model. Similar to (28), the projection error associated to PCA, if only p principal components are used, is In order to obtain a satisfactory regression model, we aim to reduce the projection error and select the principal components accordingly.
In the following, we discuss the mean squared error of prediction, which estimates the quality of the regression model, i.e., the surrogate model. The error estimators for surrogate models are analyzed onindependent test data. It may happen that the test set is not large enough, as in our case, during the start of the greedy iteration. If c k is not much larger than m, then the model may give weak predictions due to the risk of over-fitting for the parameter groups which are not used in model training. In such situations, the error estimators are designed on the learning data where leave-one-out cross-validation is an elegant choice. As we have a data set with c k observations, the leave-one-out cross-validation approach trains the model on c k − 1 data points, and test data contains one observation. We repeat this process c k times and compute the mean square error. The major advantage of this approach is that it has far less bias as if we use the entire data for training. However, the computational cost of one full leave-one-out cross-validation is high, and it may produce inconsistent and varying results [58]. Additionally, it tends to have high variance and generates very different estimates if repeated estimates with different initial samples of data from the same distribution. K-fold cross-validation is an alternative advocated to circumvent these drawbacks [20]. This technique involves randomly dividing the data set into K subsets known as folds of approximately equal size. One fold is kept for testing, and the model is trained on the remaining K − 1 folds. The process repeats for K iterations each time with K − 1 folds as the training set and the Kth fold as a test set. As we repeat the process K times, we get K errors per each fold, which are then used to compute the mean square error of prediction. The Kfold cross-validation is computationally more efficient than the leave-one-out cross-validation as we repeat the process only for K iterations instead of c k iterations.
Let L = [P k ,ε] ∈ R c k ×(m+1) be a learning data matrix consisting of c k parameter groups and error estimator values at kth iteration of the adaptive greedy algorithm. HereP k = [ρ 1 , · · · , ρ c k ] ∈ R c k ×m andε = [ε 1 , · · · , ε c k ] ∈ R c k ×1 . Let the dataset L be partitioned into K groups or folds L 1 , · · · , L K such that L i ∩ L j = ∅ for any i = j. Without loss of generality, we suppose that c k is a multiple of K where c k denotes the size of dataset L. The size of each fold is then c k /K. We train our surrogate modelε on L \ L K , i.e., all data available in the matrix L except in L K , while the data set L K is used as a test set. After K iterations, we take the sum of errors and calculate its mean to get the MSEP K−CV . The K-fold cross-validation error is then given as [40] where the inner sum is calculated over the Kth fold. The bias of . In K-fold cross-validation, the surrogate model is trained on L \ L K data instead of all L. It can be expected that this surrogate model may be inferior as compared to the model trained on all data L. Thus, it is suggested to use the adjusted K-fold cross-validation to avoid this problem. The adjustment is [40] where MSEP app. is called the apparent MSEP or mean squared error of calibration (MSEC), which uses the whole learning data L as a test set. The second term of (35) is calculated as follows. In contrast to (34), instead of testing the model trained on the excluded data set L k , we test it on the training data set L \ L k itself. After K iterations, we compute the mean by dividing the sum of K errors by c k , which is nothing but the second term of (35). The adjustment is then the difference between the apparent error and the error calculated by training and testing the model on the L \ L K data. In short, this adjustment considers the errors associated with the K folds, which are not used for the training purpose during the K-fold cross-validation. The apparent MSEP is given as Thus, the adjusted K-fold cross validation error estimate is We can now define the total sampling error associated with the PCR technique as a summation of the projection error and the adjusted K-fold cross-validation error Algorithm 6 presents a methodology to obtain a suitable surrogate model that minimizes the sampling error. It starts by computing a projection error ǫ PCA associated with the PCA. The surrogate modelε is then constructed using j = p principal components. To assess the quality of the surrogate model, we determine the cross-validation error. The data set L is divided into K folds. We train the designed surrogate model on all data available in the matrix L except in L k . Furthermore, the algorithm evaluates the surrogate model on the omitted data L k , then computes errors, and subsequently determines the K-fold cross-validation error. MSEP K−CV is then adjusted and added to ǫ PCA to get the sampling error ǫ samp . If the sampling error is less than the user-defined tolerance e samp tol , then the algorithm terminates. Finally, we use the surrogate model built using j = p principal components as explained in [7]. Build a surrogate modelε using p = j principal components.
Divide L into K random folds. Compute squared errors

Sensitivity Analysis
In the previous subsections, we addressed the errors ǫ h , ǫ POD , ǫ RM , ǫ samp associated with each stage of the model order reduction approach. This section presents a sensitivity analysis for these errors that provides quantitative information regarding the relative contribution of each error to the model output. Sensitivity analysis can be helpful in various situations, including forecasting as well as identifying where improvements or adjustments are needed in a process [14]. Such information helps to improve the model output by selecting the appropriate grid size h, the reduced dimension d, and the optimal sampling parameters c. The sensitivity indices rank the numerical errors according to their contribution to the total error. We can use this ranking to allocate the computational resource efficiently. More computational resources could be provided to the most sensitive parameters [35].
In general, there are two types of sensitivity analysis, i.e., local and global [10,14]. The local sensitivity analysis evaluates changes in the output with respect to the variations in a single factor by keeping other factors constant [30]. Such sensitivity is often evaluated through gradients or partial derivatives of the output functions at these factors. The main limitation of the local sensitivity analysis is that it evaluates the factors one at a time and does not consider the interactions and simultaneous changes in all factors [51]. Global sensitivity analysis overcomes this drawback by quantifying the importance of all factors and their interactions to the model output. It provides an overall influence of factors on the output in contrast to a local view of partial derivatives as in local sensitivity analysis.
It is easy and computationally feasible to compute the local sensitivity of deterministic errors such as the discretization error by computing the norm of the first-order factor in a Taylor series expansion or the partial derivative using finite differences [14]. For stochastic errors like the sampling error, it is more appropriate to compute the change in the output with respect to the variations in inputs [35]. However, this creates a difficulty in comparing the relative contributions of the deterministic errors versus the stochastic errors. Thus, it is necessary to propose an approach that can be applied to both types of errors.
There are different local and global sensitivity analysis methods, such as the one-at-atime (OAT) approach [10], the automatic differentiation technique (AD) [15], variance-based sensitivity analysis (the Sobol method) [37,51], and Fourier amplitude sensitivity analysis (FAST) [39]. The local sensitivity analysis approaches like OAT and AD calculate the firstorder partial derivatives of the outputs with respect to small changes in the inputs. In AD, a computationally efficient computer code automatically evaluates the partial derivatives and can be applied without having detailed knowledge of the algorithm implemented in the model. However, the technique may be limited to specific computer languages and requires specific libraries, and the accuracy for sensitivity results depends on the numerical method used in the AD software. Both the Sobol method and the FAST method are variance-based sensitivity analysis techniques capable of computing first-order (local) and global sensitivity indices. The main difference between these two approaches is the underlying algorithm in multidimensional integration of sensitivity indices. The FAST approach uses a pattern search technique based on a sinusoidal function, while the Sobol approach uses a Monte Carlo technique to determine the sensitivities. The disadvantage of the FAST approach is that it demands specific sampling and may perform poorly for discrete inputs and models with discontinuities. Among the discussed techniques, the Sobol method is one of the most powerful techniques being model-independent and robust, and it considers interactions between different input factors. However, it may become computationally intensive if there are many input factors as it uses Monte Carlo techniques. In our case, the number of input factors (numerical errors) is small, which eases the computational burden. Moreover, [35] showed that the variance-based sensitivity analysis approach could be applied considering both deterministic and stochastic numerical errors. Thus, in this work, we use the Sobol method to perform the sensitivity analysis.
Consider a model output Y depending on some input factors (X 1 , . . . , X p ) Y = f (X 1 , . . . , X p ). Now, the variance based first order effect of X i is [37] Var where X i is the ith factor and X ∼i is the p − 1 dimensional space composed of all factors except X i . The inner expectation operator means that the mean of Y is taken over all possible values of X ∼i while the factor X i is kept fixed. Var X i (·) shows that the variance is taken over all possible values of X i . Using the law of total variance, we can write [41] Var by normalizing where the first order (local) sensitivity index for the factor X i is given by the first term of (41) Here S i is a normalized index as Var X i E X ∼i (Y |X i ) varies from zero to Var(Y ) and measures the first order effect of X i on the output Y . E X i Var X ∼i (Y |X i ) is the residual [51]. It is necessary to note that the low value of the first-order sensitivity index does not mean that the factor is not important. It might make a more significant contribution to the output by interactions with other factors. Thus, the total sensitivity index is necessary, which considers the main effect of the factor and the effect of its interactions with other factors on the output. In the sensitivity analysis framework, functional decomposition of the total variance Var(Y ) is referred to as functional analysis of variance (ANOVA) and given as [54] Var where Var i1,...,ip−1 .
The first-order sensitivity index (42) can then be obtained from the first p terms of the decomposition (43) as The other terms of the decomposition can be interpreted as the higher-order sensitivity indices. The second-order index S ij represents the effect of interactions of factors X i and X j on the output and so on up to order p. Therefore, for p factors, we have 2 p − 1 sensitivity indexes. Based on the decomposition (43), [25,52] introduced the so-called total indices or total effects as where #i denotes all the indices associated with the factor X i . The total sensitivity is derived using the expressions in (43) as [51] S Using the law of total variance, similar to (40), we can write the total variance of Y by exchanging X i and X ∼i as Dividing (44) by Var(Y ) we obtain Var(Y ) .

Calculation of Sobol Indices
The methodology to determine the local and global Sobol indices with a single simulation is presented in [51]. Consider two independent sampling matricesX andX to compute the values of the output Y corresponding to different input factors X 1 , ..., X p .
where n is the sample size and the p is the number of input factors. Then each row of the matricesX andX is then a one parameter sample set. We introduce a new matrixX iX whose columns are from the matrixX expect the ith column, which is from the matrixX.
Based on these matricesX,X, andX iX , we determine the local sensitivity and the global sensitivity indices as follows. We obtain a variance of Y as in [41] Var Applying (45) to From [50] we then have and Using (46) and (47), we define the local sensitivity index as [51] The global sensitivity index S T i is derived in a similar way, see [32,51] as where f (X), f (X), and f (X iX ) are the model output results using the parameter sample values from matricesX,X, andX iX , respectively. The independent distribution of each numerical error is determined as follows. The discretization error ǫ h depends on the grid size h used to construct the full model. We sample different grid sizes with an appropriate grid refinement ratio and construct different full models. We further calculate the discretization error associated with each model and use generated samples to construct the distribution for ǫ h .
Similar to the discretization error, the projection error ǫ POD and the reduced model error ǫ RM can be regulated by varying the reduced dimension d. We have to solve the full model for some training parameter groups to construct the snapshot matrix in order to generate the reduced basis. Before obtaining the distributions of the projection error and the reduced model error, the full model dimension and the number of training parameters are fixed. We construct the reduced models with different reduced dimensions and obtain the corresponding errors ǫ POD , and ǫ RM . These error samples are then used to construct the desired distributions.
The training parameters used to construct the surrogate model govern its quality. We randomly sample the parameter groups that initiate the adaptive algorithm (parameter set P 0 ) as well as vary the cardinality of setP. The first few samples of the surrogate model error are obtained by randomly sampling the parameter setP 0 while keeping the cardinality ofP constant. Later on, the cardinality ofP is varied whileP 0 is fixed by specifying the particular seed number.
The obtained samples are then randomly separated into two matricesX, andX. Each column of the matricesX orX represents the distribution of one numerical error, while each row describes a set of model factors, i.e., in this case, the numerical error {ǫ h , ǫ POD , ǫ RM , ǫ samp }. The total error (11) presented in Section 3 is then used as a model output to perform the sensitivity analysis based on the generated sample data.

Model Hierarchy: Stopping Criteria
We now use the described techniques to calculate the financial risk associated with an invested product. Underlying risk factors such as interest rates have a direct influence on the invested asset. Interest rate risk arises when unanticipated developments shift the yield curve or change its shape. Thus, to calculate the financial risk, the instrument is evaluated for several thousand different simulated yield curves. One can consider these simulated yield curves as different scenarios. The financial instruments are evaluated via the dynamics of short-rate models M, based on the convection-diffusion-reaction partial differential equations (2). The choice of the short-rate model depends on the underlying financial instrument. The model hierarchy simplifies the process of obtaining a reduced model and is depicted in Fig. 2. It starts by discretizing the partial differential equation using either the finite difference method or the finite element method. The discretized model is the full model M h . To obtain a reduced model M d h , we compute a reduced basis via the POD approach along with the classical or adaptive greedy sampling. Finally, the full model is projected onto the reduced basis to get the reduced model. The above-described steps finally lead to the desired reduced model. However, it may happen that the user does not have enough knowledge or background about these methods to obtain the optimal reduced model. Thus, the user prefers an automatic approach that produces a reduced model whose error is below a user-defined tolerance. We describe the methodology to obtain a reduced model as a black box for which the user only requires to provide inputs and gets the desired output. The black box method automatically selects the partial differential equation according to the underlying financial instrument, generates the discretized model, obtains the training parameters, generates the desired reduced basis, and produces a reduced model with the desired user-defined tolerance. Yield curves and a financial instrument are fed as inputs to the black box, which outputs the scenario results. The procedure is illustrated in Fig. 2. This approach also yields a well defined stopping criteria as follows.
Consider a parameter space P ∈ R N ×m (of scenarios N ) for which the underlying financial instrument is to be solved using an appropriate short-rate model M. Let ℓ ∈ P ⊂ N be the number of training parameters selected using a sampling algorithm. Now, these ℓ parameters are used to obtain the reduced model M d h with the total error ǫ T ≤ ǫ h + ǫ d + ǫ samp .
Remark 3.1. The stopping criteria to valuate an underlying financial instrument with N scenarios states that when ǫ T ≤ e tol , where e tol is the user-defined tolerance, only ℓ scenarios need to be evaluated using the full model M h , while the remaining N − ℓ scenarios can be evaluated using the reduced model M d h .
We employ the ℓ scenarios with the full model M h to obtain the reduced basis and the

Numerical Example
Steepeners are the financial instruments where the cashflows depend on differences in the evolution of two different reference rates with different tenors in one currency [1]. We consider a callable/puttable steepener instrument whose coupons depend on the difference between two constant maturity swap (CMS) rates, i.e., CMS10 -CMS2 [11]. Table 2 shows the properties of a puttable steepener. A puttable steepener gives the right to the investor to claim an early redemption on certain dates (typically coupon dates) for a price fixed (typically the nominal value) in the term sheet of the bond. These early exercise dates are before the maturity of the bond. In short, a put is an option that gives the right to a buyer to sell the underlying asset at an agreed price at some time. The selling price is known as a strike price. The payoff of an option is its value at the time of its exercise. Consider K as the strike price of a put option and an underlying instrument with a value V and maturity T . The payoff of the put option is Table 2 shows that the interest rate is capped at 3.0% p.a. and floored at 0.0% p.a. The cap is nothing but an upper limit on the interest rate for a floating interest rate instrument, while the floor sets a lower limit on the interest rate. The coupon rate from the fourth year until maturity is coupon rate =min(coupon cap, max(coupon floor, CMS10 − CMS2)), =min(3.0%, max(0.0%, coupon floor, CMS10 − CMS2)) We solve the puttable steepener example using the two-factor Hull-White model. The partial differential equation is discretized using the finite element method to obtain the full model, which is further approximated to achieve the reduced model. The model parameters are σ 1 = 0.0035, σ 2 = 0.008, γ = 0.65, α = 0.75, b = 0.04. We have computed the model parameter θ ℓ as explained in [38]. The yield curve simulation is the first step to compute the model parameters. Based on the yield curve simulation procedure described in [7], we have performed the bootstrapping process for the recommended holding period and the intermediate holding period, i.e., for ten years and five years, respectively. The collected historical data has 21 tenor points and 1306 observation periods as follows  Figure 3 shows the plots of simulated yield curves obtained by bootstrapping for five years and ten years in the future. The calibration generates a parameter space P of order 10 000×m for which we now solve the two-factor Hull-White model. For calculations, we have used a structured, two-dimensional triangular mesh grid that generates a full model of order M . We have used fully implicit time-stepping so that discretizations in time and space can be chosen independently. The resulting system of linear equations is To obtain a solution, we solve this system starting at t = T with an appropriate terminal condition V (T ) backward in time at each time step n. It is important to note that we have selected the time step ∆t = 20 days using Algorithm 1, and ensuring that all key dates are achievable, i.e., coupon dates, put dates, or valuation dates. The tolerance e t tol of 10 −3 is used to obtain this best suitable time step ∆t. We determine the value of the steepener by solving the linear system given in (49) as follows. If the coupon date is reached, then we update the value of the steepener by adding the necessary coupon given by (48) [38] where V n−1 is the value calculated by solving the full model, and C F is either an already known cashflow or a cashflow generated by the coupon if the coupon date is achieved at t n−1 or zero. Similarly, if the put date is reached, then we have to update the value at that time point by the following max function where the put price is the strike price defined by the put option contract. In this paper, we have considered Put Price = 1. We have considered put dates annually starting one year after the issuing of the bond until one year before maturity. The quality of the full model is assured by calculating the discretization error associated with the two-factor Hull-White model, as explained in Subsection 3.1. Figure 4 shows the discretization error along with the grid convergence index values. We can see that the error decreases with increasing order of the full model. A discretization error computation requires three different solutions of the full model for three different grid sizes. The error estimator is acceptable if these solutions are in the asymptotic range of convergence, i.e., the observed order of accuracy matches the formal order of accuracy, and the grid convergence index is low. The formal order of accuracy depends on the order of the polynomial used to define the elements in the finite element method [29]. The pth order polynomial gives the p + 1 order of convergence [4]. We have constructed our full model using the 1st order linear triangular elements. Thus, the formal order of accuracy p f = 2. The observed order of accuracyp is in the range of 1.8959 to 2.5616 with an average of 2.072. For the full model with M = 1600, p = 1.9213, which matches the formal order of accuracy within 10%. Thus, we can say that the solutions obtained using the full model are in the asymptotic range. The safety factor used to calculate the grid convergence index (26) is then 1.5. The full model size is selected using Algorithm 5, which is terminated when the discretization error is less than the userdefined tolerance e h tol = 5 × 10 −4 . We have designed a full model of order M = 1600 as the discretization error ǫ h = 2.182 × 10 −4 falls below 5 × 10 −4 . Also, one can notice that the grid convergence index is less than 0.05%.

Model Order Reduction
For the model order reduction approach, we have obtained the training parameter set to compute the optimal reduced basis by using either the classical greedy sampling or the adaptive greedy sampling. In the following, we present the results and error analysis of these algorithms. For the classical greedy sampling, we have used c = 40 randomly selected parameter groups to train the algorithm. The allowed maximum iterations I max = 10, i.e., at most 10 parameter groups are selected which maximize the residual error. The algorithm then updates the snapshot matrix and generates a new reduced basis.  Figure 5 shows the evolution of maximum and average residuals with each iteration of the classical greedy sampling algorithm. It is observed that the residual error decreases with each proceeding iteration. We can truncate the greedy iteration after 4th or 5th iteration as the error has dropped below 10 −4 .
The dimension of the reduced basis d is chosen using Algorithm 5. The projection error is plotted in Fig. 6. We observe that the POD projection error decreases monotonically with increasing d, which is determined in terms of POD eigenvalues. The graph of the projection error shows that we have succeeded in determining an optimal reduced basis, which is later used to generate the reduced model. Figure 6 also shows the reduced model error ǫ RM is trained on the randomly generated parameter setP of cardinality 40. We noticed that the random sampling often neglects some of the vital parameter groups within the parameter space. The reduced model designed using the classical greedy sampling approach shows this drawback as illustrated in Fig. 7. The figure shows the reduced model error for two different parameter groups (ρ 342 , ρ 5287 ). In both cases, the relative error does not drop below the user-defined tolerance e d tol = 5 × 10 −4 , and Algorithm 5 does not converge. This prompted us to implement an adaptive greedy approach that relies on the surrogate model to locate the training parameters. The used surrogate model is based on the principal component  regression technique [7]. The adaptive greedy iteration is initiated by computing the error estimator for randomly selected parameter setP 0 of cardinality c 0 = 10. At each surrogate model iteration, we have chosen c k = 10 parameter groups. This process repeat itself until we get the parameter setP of cardinality 40. Figure 8 shows monotonically decreasing singular values of the data matrixP k composed of c k parameter groups at kth iteration. We noticed that the first three principal components comprise almost 99% of energy, and within the first five components, energy level round up to 99.99%. Thus, the surrogate model has only five principal components. Figure 8 demonstrates that the predicted values obtained via the surrogate model match the actual values within an acceptable range. The quality of the designed surrogate model is further tested using a mean squared error of prediction. We have used a K-fold cross-validation technique to determine the MSEP as described in Subsection 3.3. The learning data matrix composed of k × c k parameter groups and error estimators at the kth iteration is divided into K = 4 folds. We obtain the mean square error for each fold and sum it to calculate the MSEP of the cross-validation. Table  2 shows the calculations and MSEP for the designed surrogate model. The sampling error ǫ samp is the sum of the adjusted K-fold cross-validation error and the projection error, i.e., ǫ samp = MSEP K−CVadj + ǫ PCA = 2.3750 × 10 −4 . As ǫ samp < e samp tol = 5 × 10 −4 , we can say that the designed surrogate model is satisfactory. We have used this surrogate model to construct the parameter setP. Finally, the optimal parameter group ρ I is selected from the parameter setP that maximizes the residual error. The reduced basis is obtained by computing the truncated SVD of this snapshot matrix, as shown in Algorithm 5. Figure 9 shows the evolution of the maximum and average residual error with each adaptive greedy iteration. It illustrates that the residual error decreases with increasing iteration and falls below 10 −4 after the fourth iteration, which is enough to truncate the greedy iteration. The decrease of the residual error with each incrementing iteration also shows that the algorithm succeeded in efficiently locating the optimal parameter group. The reduced basis dimension d is selected such that for the generated reduced model, ǫ POD and ǫ RM are below the user defined tolerance. Figure 10 shows the monotonically decreasing projection error ǫ POD associated with the POD. The reduced model error ǫ RM plotted in Fig. 10 shows that the error decreases with increasing the reduced dimension d. We have constructed the final reduced model using d = 10, as the model order reduction error is less than the user-defined tolerance e d tol , i.e., ǫ d = ǫ POD + ǫ RM = 4.042 × 10 −6 + 1.882 × 10 −4 = 1.922 × 10 −4 < 5 × 10 −4 . This means that the reduced model is an excellent approximation of the full model. We have performed 10 full model evaluations during the adaptive greedy sampling, and the remaining scenarios are solved using this reduced model, as explained in Subsection 3.5. We have implemented the adaptive greedy sampling approach to overcome the drawbacks of the classical greedy sampling method, as evident from Fig. 11. The relative error between the full model and the reduced model obtained using the adaptive greedy sampling approach is smaller than 10 −4 with the reduced dimension d = 10 for both parameter groups ρ 342 , ρ 5287 . The error ǫ RM keeps on decreasing with increasing d, unlike in the classical greedy sampling approach.
The convergence of the classical greedy sampling approach is monitored using the residual error. However, we have used an error modelǭ RM that uses the reduced model errors and the residual errors to monitor the convergence of the adaptive greedy sampling. Figure 12 illustrates the designed error model based on the available error set E p for 4 different greedy iterations. The error plot demonstrates a strong correlation between the relative error and the residual error. The results indicate that a consideration of the linear error model is adequate to capture the overall behavior of the exact error as a function of the residual error.

Random Sampling
We have used random sampling up to a certain extent in both classical and adaptive greedy approaches. Random sampling is used to construct the parameter setP in the classical greedy approach, while in the case of the adaptive greedy approach, it is used to initiate the algorithm. In both cases, the random sampling somewhat affects the results. Figure 13 shows the maximum residual error plot with each proceeding iteration of the greedy algorithms with error bars. These results are obtained by performing classical and adaptive procedures for 30 different sets of random numbers. It is observed that the adaptive greedy approach is less sensitive to random sampling than its classical counterpart. On average, the adaptive greedy approach yields a better reduced model as its residual errors are smaller than the classical approach.

Sensitivity Analysis
We have obtained 30 different data points for each numerical error to generate respective distributions by varying certain key parameters. We have sampled 30 different grid sizes to compute the discretization error distribution, where the full model dimension is being varied from 100 to 2500. To compute the sampling error distribution, we have randomly generated 15 differentP 0 by keeping the cardinality of n(P) = 50. Additionally, 15 data points for the sampling error are calculated by varying n(P) between 30 and 100, while keepingP 0 fixed using a particular seed number. The distributions for the projection error and the reduced model error are obtained by varying the reduced dimension between 5 and 60.
The local indices show that the final output is most sensitive to the discretization error followed by the sampling error, while it is least sensitive to the reduced model error and the projection error, respectively. We can mitigate the discretization error by increasing the full model dimension. Nevertheless, sensitivity indices prevail that the change in grid size considerably affects the final output. We overcome this problem by increasing the full model dimension until the further change in grid sizes does affect the solutions significantly. We have constructed the full model of order M = 1600, where the discretization error is of order 10 −4 . The global sensitivity analysis shows that the sampling error is the most sensitive numerical error. The major input factors contributing to the sampling error are the randomly sampled parameter setP 0 and the cardinality of setP. As it is not feasible to search the entire parameter space as the parameter setP, the surrogate model only locates a certain number of parameter groups to constructP and neglects the rest. This omission of parameter groups contributes to the sampling error together with the random sampling ofP 0 used to initiate the surrogate model loop, as we have already demonstrated in Figure 13. This random sampling contributes to the uncertainty and makes the sampling error more sensitive. To overcome this problem, we suggest to increase the cardinality ofP used to locate the most optimal parameter group if necessary. We have kept the cardinality ofP as 40, as the sampling error is of order 10 −4 .
Although increasing the cardinality of the set will increase the computational burden, it generates the best suitable reduced basis that allows constructing a fairly low dimensional reduced model. This compensates for the extra computational burden necessary during the adaptive sampling. We noticed that the projection error is the least sensitive to the final output considering both local (individual contribution) and global (interactions) indices. This tells us that the singular value decay of the snapshot matrix is efficient to obtain the desired reduced basis.

Computational Cost
The three major error contributors that lead to the total error are the discretization error, the model order reduction error, and the sampling error. Figure 14 presents the model hierarchy  along with the error associated with each stage. The total error is a sum of all significant errors and is given by (11) ǫ T ≈ ǫ h + ǫ RM + ǫ POD + ǫ samp. , ǫ T = 2.182 × 10 −4 + 1.882 × 10 −4 + 4.042 × 10 −6 + 2.3750 × 10 −4 = 6.479 × 10 −4 .
As ǫ T < 10 −3 fulfills the stopping criteria prescribed in Remark 3.1. Thus, we can say that ℓ = 10 full model valuations are enough to obtain a fairly accurate reduced model, and the remaining scenarios (parameter groups) can be solved using the generated reduced model.
Calculation of the computation time is a crucial aspect enlightening the importance of the model order reduction over the full model approach. In the classical greedy sampling approach, the algorithm solves c reduced models and one full model at each greedy iteration. The algorithm then updates the snapshot matrix using the full model solution and computes its truncated SVD. Let t RM be the time required to solve one reduced model, t FM be the computational time required for one full model, and let t SVD be the time required to obtain a truncated SVD of the snapshot matrix. The total computational time T CG Q required to obtain the reduced basis after i iterations of the classical greedy sampling approach can then be given as The adaptive greedy approach at each iteration conducts the following steps to obtain a reduced basis. It is initiated by solving c 0 reduced models to obtain the residual errors {ε j } c 0 j=1 . These error estimator values are used to construct a surrogate model that locates c k parameter groups. The process is repeated until the total cardinality ofP reaches c, i.e., for k iterations. The algorithm then solves the full model for the optimal parameter group, updates the snapshot matrix, and computes a truncated singular value decomposition of the updated snapshot matrix. Finally, an error modelǭ RM is built to monitor the convergence of the greedy algorithm. Thus, the total computational time can be given as T AG Q ≈ c 0 × t RM + k(c k × t RM + t SM + t ev SM ) + t FM + t SVD + 2t af t,bef where t SM and t ev SM denote the computational times required to build and valuate a surrogate model for the entire parameter space, respectively where t EM is the time required to build an error model. The term 2t af t,bef RM shows the computational time needed to solve the reduced model after and before updating the reduced basis. We have used a truncated SVD based on randomized algorithms to obtain a reduced basis. Table 4 shows the computation time comparison between the randomized truncated and full SVD. For our problem, the randomized SVD is at least 12 times faster than the basic SVD. We noticed that the evaluation of the reduced model is at least 8 − 10 times faster than the full model. However, there is a slight increase in the total time due to the additional reduction time. Table 6 also infers that the evaluation time increases with an increasing reduced dimension d. The reduced model obtained based on the classical greedy sampling approach is at least 7 − 9 times faster than the full model. The computational time for the adaptive greedy based reduced model is a bit higher due to the time invested in computing surrogate and error models. The reduced  model obtained using the adaptive greedy sampling approach is 6 − 8 times faster than the full model. The computation time presented in the table considers that the greedy algorithms run for the maximum number of iterations I max . However, in practice, we can truncate the algorithms after the 4th or 5th iteration and achieve even more speedup.

Steepener Results
We solve the steepener instrument with the reduced model obtained using the adaptive greedy approach for 10 000 different parameter groups. To design a key information document, we need the values of the steepener at different spot rates. The spot rate is nothing but the yield rate at the first tenor point from the simulated yield curve. For 10 000 simulated yield curves, we procure 10 000 different spot rates and the corresponding values for the instrument. The regulations also demand to find a VaR equivalent volatility (VEV) to determine a market risk indicator at the recommended holding period. The VaR shall be the instrument's price at a confidence level of 97.5% at the end of the recommended holding period discounted to the present date using the risk-free discount factor (DF). Thus, the VaR in price space is VaR × DF, and the VEV is VEV = (3.842 − 2 × ln(VaR price space )) − 1.96 √ T .
For the steepener instrument, the calculated VEV is 1.3285% which gives the market risk Figure 15: Distribution of 10 000 results after five years (left) and ten years (right).
indicator as 2. See [7] for a detailed discussion on the VEV and market risk indicators. It is required to include three different performance scenarios in the KID: (i) favorable scenario, (ii) moderate scenario, (iii) unfavorable scenario, which are the values at 90th percentile, 50th percentile, 10th percentile of 10 000 values, respectively. The steepener performance scenarios obtained using the reduced model and the commercially available software UnRisk for the comparison are presented in the table 7 after five and ten years. Figure 15 shows the distribution of fair values of the instrument, plus the coupons in the respective path obtained so far after five years and ten years in the future. The histogram gives an idea about the frequency distribution and shape of a set of 10000 values for the PRIIP.

Conclusion
This paper presents a detailed error analysis of the model order reduction framework developed in [7] for the financial risk analysis. Three major sources of errors considered are the discretization error, the sampling error, and the model order reduction error. We have developed error estimations and their quantification methods for each error. Based on these error estimators, we have designed algorithms to select optimal grid size, suitable training parameters, and the appropriate reduced dimension. The algorithms are designed such that the model order reduction framework can be used as a black box for which the user only requires to provide inputs and get the desired output. Another major contribution of this paper is the implementation of a sensitivity analysis to rank the contribution of each numerical error. We have developed a methodology to address both deterministic (e.g., discretization error) and stochastic errors (sampling error).
The Sobol variance-based sensitivity analysis is used, which computes local as well as global indices. The global sensitivity analysis revealed that the sampling error is most sensitive to the final output. Such information advised us to allocate more computational resources to the sampling algorithm. This ultimately generated a suitable reduced basis as demonstrated by the obtained results.
We have tested the designed algorithms for a numerical example of a puttable steepener under the two-factor Hull-White model. The reduced model is obtained such that the total error ǫ T = 6.479 × 10 −4 is less than the user-defined tolerance e tol = 10 −3 . As ǫ T < e tol , we can conclude that the full model valuations for the adaptively selected ℓ = 10 parameter groups are enough, and the rest of the parameter groups can be solved inexpensively using the reduced model of dimension d = 10. The reduced model approach to solve the entire parameter space was at least 8−10 times faster than the full model approach. We also noticed that the randomized singular value decomposition provides an excellent speedup compared to the basic singular value decomposition.
We conclude that the new model order reduction approach shows many potential applications in the historical or Monte Carlo value at risk calculations.