Yield curve interpolation using constraint kriging techniques
Due to the lack of reliable market information, building financial term-structures may be associated with a significant degree of uncertainty. In this paper, we propose a new term-structure interpolation method that extends classical spline techniques by additionally allowing for quantification of uncertainty. The proposed method is based on a generalization of kriging models with linear equality constraints (market-fit conditions) and shape-preserving conditions such as monotonicity or positivity (no-arbitrage conditions). %and where interpolation constraints are given in terms of linear systems. %\hassansuppr{The curve we propose to use for financial application corresponds to the most likely path (the mode) of the underlying conditional Gaussian process. We define the most likely curve and show how to build confidence bands. The Gaussian process covariance hyper-parameters under the construction constraints are estimated using cross-validation techniques. Based on observed market quotes at different dates, we demonstrate the efficiency of the method by building curves together with confidence intervals for term-structures of OIS discount rates, of zero-coupon swaps rates and of CDS implied default probabilities. We also show how to construct interest-rate surfaces or default probability surfaces by considering time (quotation dates) as an additional dimension.
Kriging_of_Financial_Term_Structure_06_avril_2016.tex
\documentclass[12pt,a4paper]{article}
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{graphicx}
%did
\usepackage[usenames,dvipsnames]{xcolor}
%\usepackage{xcolor}
\usepackage{subfig} 
\usepackage{lmodern}
\usepackage[left=2cm,right=2cm,top=2cm,bottom=4cm]{geometry}
\usepackage{float} % Pour fixer la position d'une figure
\usepackage{dsfont} % pour indicatrice
\usepackage[pdftex,colorlinks,citecolor=blue, urlcolor=blue]{hyperref}
%\RequirePackage[colorlinks,citecolor=blue, urlcolor=blue]{hyperref}
\usepackage{natbib}

%\usepackage{placeins}%pour \FloatBarrier

\newcommand{\cov}{\mathrm{Cov}}
\newcommand{\var}{\mathrm{Var}}
\newcommand{\indic}[1]{\mathds{1}_{#1}}
\newcommand{\Nx}{N_x}
\newcommand{\Nt}{N_t}
\newcommand{\Dom}{D}
\newcommand{\ta}{\tau}
\newcommand{\suchthat}{\;\ifnum\currentgrouptype=16 \middle\fi|\;}

%\def \proof{{\noindent \bf Proof. }}

\usepackage[normalem]{ulem}
\newcommand{\hassan}[1]{{\color{blue}#1}}
\newcommand{\hassansuppr}[1]{\hassan{\xout{#1}}}

%\newcommand{\ares}[1]{{\color{red}#1}}
%\newcommand{\did}[1]{{\color{BlueViolet}#1}}
%\newcommand{\didsuppr}[1]{\did{\xout{#1}}}

\newtheorem{theorem}{Theorem}[section]
\newtheorem{lem}[theorem]{Lemma}
\newtheorem{pro}[theorem]{Proposition}
\newtheorem{cor}[theorem]{Corollary}
\newtheorem{conj}[theorem]{Conjecture}
\newtheorem{rem}[theorem]{Remark}
\newtheorem{com}[theorem]{Comments}
\newtheorem{ex}[theorem]{Example}
\newtheorem{defi}[theorem]{Definition}
\newtheorem{hyp}[theorem]{Assumption}

\def \CC{{\mathbb C}}
\def \DD{{\mathbb D}}
%\def \\E\E{{\mathbb \E}}
\def \FF{{\mathbb F}}
\def \GG{{\mathbb G}}
\def \HH{{\mathbb H}}
\def \LL{{\mathbb L}}
\def \MM{{\mathbb M}}
\def \NN{{\mathbb N}}
\def \PP{{\mathbb P}}
\def \R{{\mathbb R}}
\def \ZZ{{\mathbb Z}}

\def \Cxi{{\mathcal C}}
\def \Ixi{{\mathcal I}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\annu#1{_{% 
  \vbox{\hrule height .2pt 
    \kern 1pt 
    \hbox{$\scriptstyle {#1}\kern 1pt$}% 
  }\kern-.05pt 
  \vrule width .2pt 
}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------------------ gestion des graphiques
\usepackage{ifpdf}
\ifpdf
% useless with Elsarticle package:
 % \usepackage[pdftex]{graphicx}
 % \DeclareGraphicsExtensions{.pdf,.jpg,.png}
%non utilisé:
  %\pdfcompresslevel=9
  %\usepackage[pdftex]{hyperref} %hyperref permet les liens hypertextes dans le pdf
\else
  \usepackage{graphicx}
  \DeclareGraphicsExtensions{.eps,.jpg,.png}
\fi
\usepackage{epstopdf}
%------------------------------------ fin de gestion des graphiques

\title{Yield curve interpolation using constraint kriging techniques}

\begin{document}
\begin{center}
\Large{\textbf{Kriging of financial term-structures}}
\end{center}
%\tableofcontents

\vspace{0.4cm}

\begin{center}
 Areski Cousin$^{\dagger}$\footnote{areski.cousin@univ-lyon1.fr}, Hassan Maatouk$^{\ddagger }$\footnote{hassan.maatouk@mines-stetienne.fr} and Didier Rulli\`ere$^{\dagger}$\footnote{didier.rulliere@univ-lyon1.fr}
\end{center}
 
\vskip0.2cm

\begin{center} 
($\dagger$) 
{\small Universit\'{e} de Lyon, Universit\'{e} Claude Bernard Lyon 1, ISFA, Laboratoire SAF  EA2429, 50 av. Tony Garnier, 69366 Lyon, France\\}
%{\small Universit\'{e} de Lyon, Universit\'{e} Lyon 1, ISFA, LSAF, 50 avenue Tony Garnier, 69366 Lyon, France\\}


($\ddagger$) {\small Mines de Saint-Etienne, 158 Cours Fauriel, 42023 Saint-\'Etienne, France\\}
\end{center}
%
%\vskip2cm


\begin{abstract}
Due to the lack of reliable market information, building financial term-structures may be associated with a significant degree of uncertainty.
In this paper, we propose a new term-structure interpolation method that extends classical spline techniques by additionally allowing for quantification of uncertainty. 
The proposed method is based on a  generalization of kriging models with linear equality constraints (market-fit conditions) and  shape-preserving conditions such as monotonicity or positivity (no-arbitrage conditions). %and where interpolation constraints are given in terms of linear systems. 
%\hassansuppr{The curve we propose to use for financial application corresponds to the most likely path (the mode) of the underlying conditional Gaussian process.} 
We define the most likely curve and show how to build confidence bands.
The Gaussian process covariance hyper-parameters under the construction constraints are estimated using cross-validation techniques.
Based on observed market quotes at different  dates, we demonstrate the efficiency of the method by building curves together with confidence intervals for 
term-structures of OIS discount rates, of zero-coupon swaps rates and of CDS implied default probabilities. We also show how to construct interest-rate surfaces or default probability surfaces by considering time (quotation dates) as an additional dimension.
\end{abstract}

%\paragraph{This version} October 7, 2015

\paragraph{JEL classification} C63; E43; G12


\paragraph{Keywords} Model risk; interest-rate curve; yield curve; OIS discount curve; implied default distribution; kriging; no-arbitrage constraints


\section{Introduction}

%\did{j'ai aussi supprimé des doublons dans les labels d'equation eq:market_fit_Swaps_Euribor}


%Constructing the term-structure of interest rates 
%
% %fixed-income, FX or default-risky securities 
% is at the heart of modern asset pricing and risk management. 
% 
Constructing  term-structures is at the heart of asset pricing and risk management. 
A term-structure is a curve which describes the evolution of some financial or economic quantities as a function of time horizon. Typical examples are the term-structure of risk-free interest-rates, the term-structure of bond yields or credit spreads, the term-structure of default probabilities or the term-structure of stock return implied volatilities. These curve are typically not directly observed in the market. Thus, the curve construction is based on a  benchmark set of contingent financial instruments whose values explicitly depend on some part of the curve.  In practice, market quotes of these products only provide a partial information on the term-structure since they can only be considered to be reliable  for a small set of liquid maturities.  
% it   benchmark financial instruments whose values  depend on different points of the curve. 
%For instance, this can be fixed-income products for interest-rate term-structures, credit default swaps for the term-structure of implied default probabilities, equity options for term-structure of stock return implied volatilities. 
%However, market quotes of these products are typically reliable only for a small set of liquid maturities. 
The problem is then to transform a small set of market quotes into a continuum set of values representing the evolution of the underlying quantity of interest with respect to time horizon.\\ 


On practical grounds, the curve is assumed to belong to a family of parametric functions (Nielson-Siegel functional~\cite{nelson1987parsimonious}, polynomial splines, \citet{SW:2001}) and its construction consists in finding the underlying parameters that best fits observed market quotes for all available maturities.
%Consequently, the financial industry relies on interpolation or calibration schemes to supplement market information
%such as polynomial splines in order to supplement market information on the term-structure  while extracting the quantities of interest for all required maturities.\\
In \citet{deAndresGomez2004}, the interest-rate term-structure is estimated from bid-ask spreads of underlying instruments using fuzzy regression techniques.
 \citet{Hagan2006} provide a review of different interpolation techniques for curve construction. They introduce a monotone convex method and postulate a series of quality criterion such as ability to fit market quotes, arbitrage-freeness, smoothness, locality of interpolation scheme, stability of forward rate and consistency of hedging strategies. \citet{Andersen2007} analyzes the use of hyperbolic tension splines for construction of interest-rate term structures. The underlying optimization allows the user to control the relative importance of fit precision with respect to shape preservation (smoothness of the curve, penalization of oscillations and excess convexity/concavity).  In the same vein, \citet{Chiu_el2008} shows that $L_{1}$ cubic splines minimizes the curve oscillation without sacrificing good approximation of the data.
 %Le Floch: hyperbolic tension splines to smoothly manipulate locality and shape preservation 
  %\citet{Jerassy-Etzion2010}  
\citet{Iwashita2013} makes a survey of non-local spline interpolation techniques which preserve stability of forward rates. 
%Two-dimensional extensions of 
% \citet{LeFloch2013} introduces another quality criteria related to consistency of hedging strategies. He postulates that, given a constructed term structure, the sum of sequential deltas should be close enough to the corresponding parallel delta. He observes that most of spline techniques are not able to achieve this property correctly.
 % (locality of deltas, sum of sequential deltas close enough to the corresponding parallel delta). 
%\citet{Jerassy-Etzion2010}
Other papers such as  \citet{Ametrano2009}, \citet{Chibane2009}, \citet{Kenyon2012} or \citet{Fries2013} are concerned with the adaptation of curve construction methods in a multi-curve interest-rate environment. 
%\cite{Hagan2006} and \cite{LeFloch2013}  propose the following quality criteria for a curve construction method:
%% good yield-curve construction should satisfy :
%ability to fit market quotes of instruments used to build the curve, arbitrage-freeness, smoothness, locality of the interpolation scheme, stability of forward rate, consistency of hedging strategies (locality of deltas, sum of sequential deltas close enough to the corresponding parallel delta).
%These properties are difficult to obtain simultaneously. Some of them such as smoothness and locality are even conflicting or incompatible. As a consequence, 
Note that, in terms of interpolation scheme, there is no consensus towards a particular best practice method in all circumstances. In addition, the previous  approaches does not account for the uncertainty embedded in the process of curve construction. This could be of primary importance given that the market inputs may be  unreliable or even inexistent for some  maturities.
\\

This issue is related to the study of model uncertainty and its impact on risk management. This topic has been studied since a certain time period and, following the recent financial crisis, has received a particular interest. Impact of model risk on valuation and hedging of financial derivatives have been treated by, among others,  \citet{Derman1996}, \citet{Eberlein1997}, \citet{ElKaroui1998}, \citet{Green1999}, \citet{Branger2004}, \citet{Cont2006}, \citet{Davis2007}, \citet{Henaff2010}, \citet{Morini2011}. In most papers, the question of model risk is restricted to the class of  derivative products. One of the main objective is to quantify model uncertainty, for instance  to obtain bounds for the arbitrage free value of some derivative instruments, given some information on the underlying securities, such as marginal distribution of its price at some particular time horizons.
In contrast, the question of model risk embedded in the construction of marginal distribution or term-structure function themselves has not been investigated as a main object, whatever it  may concern discount curves, zero-coupon curves, swap basis curves, bond term structures or CDS-implied survival curves. \\


%In mining, geostatistics is extensively used in the field of mineral resource and reserve valuation - the estimation of grades and other parameters from a relatively small set of borehole or other samples. Geostatistics is now widely used in many other fields. Obviously there are geological and geographical applications. However, the techniques are also used in such diverse fields as hydrology, ground water and air pollution, soil science and agriculture, forestry, epidemiology, management of wildlife and weather prediction. 

From \citet{KW:1970} and \citet{MKGL:1996}, it is well-known that spline fitting is a special case of kriging  \citep[see also][]{bay:hal-01136466, 2016arXiv160202714B}. In addition, kriging allows to account for quantification of uncertainty. Kriging has been developed in geosatistics to estimate the density of some mineral resource in the ground given a relatively small set of borehole, see \cite{krige1951statistical}, \cite{matheron1963principles}, \cite{cressie1990origins}.  Its principle relies on the determination of the conditional distribution of a spatial random field given a set of observed values of the field. The main interest of this method is that it allows to build a predictor of quantities of interest at other locations, as well as uncertainties relying on this prediction.

Kriging is now widely used in many fields like hydrology, air pollution, epidemiology, weather prediction, etc. to interpolate some quantity of interest given some known values at different locations. Despite its popularity, there are relatively few works concerning kriging in actuarial sciences or in finance. Many reference academic journal of these fields give only few or even no entry corresponding to the word ``kriging''. The method is however sometimes referred to using the terms ``Gaussian Processes'' and ``machine learning''. Some existing works using kriging methodology in actuarial sciences and finance concern for example dynamic lifetime adjustments \citep{Debon2010327}, variable annuities valuation \citep{Gan2013795, Gan2015138}, nested simulation of expected shortfall \citep{liu2010stochastic}, Vasicek model calibration \citep{sousa2012machine}, stock market linkages \citep{Asgharian20134738} or credit scoring \citep{FernandesArtes2015}. Other works using spatial techniques are~\cite{Kanevski2008} on interest rates,~\cite{Benth:2015} for energy futures prices.
Some preprints or conference papers also mention the fit of some financial models \citep{stutvoet2007fitting}, spatial insurance \citep{paulson2006spatial}, trading and hedging strategies \citep{baysal2008response}, valuation of Bermudan options \citep{Ludkovski:2015}. Kriging methods naturally rely on some assumptions on the underlying random fields, and one must carefully consider all conditions that must be satisfied before constructing a kriging model.\\

In practice, the term-structure under construction has to satisfy several type of conditions.
%Several conditions are imposed in the process of constructing a term-structure.
One of the most important condition is the compatibility of the curve with market data, i.e., if the curve is used to value a benchmark set of instruments (under a specific pricing rule), the resulting values shall be as close as possible to the observed market quotes. In many classical situations, the market-fit condition translates into a system of linear constraints which can be easily incorporated in kriging techniques.
In addition, kriging can also handle the presence of noisy observations (using the so-called \textit{nugget effect}). This may be relevant in situation where, due to the lack of liquidity, market quotes cannot be considered to be reliable. It is then possible to incorporate an additional level of uncertainty (degree of confidence) associated with market observations.
%This condition can be extended by considering noisy observations are noisy not reliable, so that the latter   has then to be considered as noisy observation.
Monotonicity constraints also appears to be important in many applications. For instance, the price of default-free zero-coupon bonds (or risk-free discount factors) is a non-increasing  function of time-to-maturities under no-arbitrage assumption. Survival functions inferred from CDS spread term-structures are $[0,1]$-valued non-increasing functions.\\

Recently, some authors have studied the integration of monotonicity constraints into Gaussian process emulators, see e.g. \cite{doi:10.1137/140976741} and \cite{kleijnen2012monotonicity}. However, these methods do not guaranty monotonicity constraints in the entire domain. The article~\cite{abrahamsen2001kriging} also deals with the introduction of constraints at some location of a Gaussian process. In \cite{maatouk:hal-01096751}, classical kriging has been improved to tackle monotonicity, positivity constraints or bounds constraints on the curve values. In the present paper, we show how ``constrained'' kriging techniques can be used to extend the classical spline interpolation approaches by additionally quantifying the uncertainty in some illiquid part of the curve. \\


%\citet{GM:2004}

%We define admissible curves as arbitrage-free term structures which are perfectly compatible with input market quotes and which admit a minimum degree of smoothness. The aim of this paper is to provide a methodology for assessing the range of these admissible curves and its impact on asset valuation and hedging.
%% analyze the diversity of these admissible curves.
%Present values of instruments involved in the curve construction process are assumed to have linear representations in terms of some elementary quantities which define the curve. Depending on the context, these elementary quantities can be discount factors, zero-coupon prices or survival probabilities. Under this assumption and using a no arbitrage argument, we present an iterative procedure which allows to compute model-free bounds for these elementary quantities. We show that admissible curves can be easily generated using mean-reverting term-structure models, which are classically  employed for valuation and hedging of financial derivatives. Contrary to an HJM approach where an arbitrary initial term-structure is given as a model input, the initial curve is here constructed as a by-product of a calibrated term-structure model. The proposed curve construction method relies on a piecewise-constant specification of the mean-reverting level. We identify under which condition the implied parameters lead to an admissible curve.\\ 

%As can be seen in the literature, there exists a variety of construction methods
%\citet{Cont2006}
%\citet{Crepey2012}

%\section{Motivations of the paper}
%One of the main purposes of a yield curve is to value and price a portfolio of instruments and financial contracts.  To value a contract we need to calculate the present value of each individual cash flow which makes up a contract.  To do this we need a discount function or a discount rate.%%In risk management, present values of financial or insurance products are computed by discounting future cash-flows. 
%In market practice, discounting is done by using the current yield curve, which gives the offered interest rate as a function of the maturity for a given type of debt contract. This enables risk managers  to price any fixed-income
%security delivering known cash flows in the future (e.g., a fixed coupon bond).
%%%This allows to consider different interest rate for different horizons and different borrowers. 
%%
%% 
% Several operational methods are used in practice to reconstruct a yield curve from observed market data such as Euribor, FRA rates and Swap rates. 
 
%When the first three conditions are satisfied then we say that the yield curve is admissible. Depending on the context we will say that an IR curves is arbitrage-free when the forward rates are non-negative forward or equivalently, the (pseudo) discount factors are nonincreasing with respect to time-to-maturities.
%The motivation behind  this research is first to explain how admissible yield curves can be constructed via mean-reverting term-structure models and then we try to analyze to which extent the range/diversity of admissible curves affect the present values of financial  assets and insurance  liabilities.


The paper is organized as follows. Section~\ref{Sec:term-structure} states the term-structure construction problem and gives specific examples of market-fit conditions and shape preserving constraints. In Section~\ref{IntModels}, we briefly recall Gaussian process modeling with interpolation conditions or with more general linear equality constraints. In Section~\ref{IntModIConst}, we present the model defined in \cite{maatouk:hal-01096751} to incorporate monotonicity constraints into a Gaussian process emulator. We then study the associated properties such as convergence to the constrained interpolation spline and the estimation of the covariance hyper-parameters.
 In Section~\ref{EmpInv}, based on real market data, we construct curves together with confidence intervals  for different financial term-structures such as OIS discount curves, zero-coupon swap curves and CDS-implied default distributions. 
% Some numerical simulations of Gaussian processes conditionally to monotonicity constraints and given real data to estimate discount factors in one and two dimensions are shown in Section~\ref{EmpInv}.
%Section \ref{sec:admissible_curves} defines what is understood as admissible term-structures, with a focus on curves constructed from market quotes of interest-rate sensitive (IR curves) and default-sensitive products (credit curves). Our framework relies on the assumption that present values of building instruments have linear representations with respect to some elementary quantities. We show that this assumption is usually satisfied for construction of bond term-structures, OIS discount curves, forward rate curves implied from market quotes of swaps versus Libor or Euribor index rates. 
%  Section \ref{sec:bounds} explains how to compute arbitrage-free bounds for OIS discount curves and CDS survival curves.
%In Section \ref{sec:construct_admissible_curve}, we develop a methodology that allows to generate admissible curves as a by-product of a mean-reverting term-structure model. By playing with some free parameters, we illustrate the range of admissible curves that can be obtained in different approaches, for OIS discount curves and CDS-implied survival curves.  Section \ref{sec:conclusion} concludes.



\section{The term-structure construction problem}
\label{Sec:term-structure}
%The market-fit and the shape-preserving condition

The main ingredients in the construction of a term-structure function is a set of  financial products whose value depends on some points of the curve. Then, observing the price of these products provides an indirect (and partial) information on the curve. The first step is then to specify the relation between the value of these products and the values of the curve at different time horizons. In this paper, we restrict ourselves to situations where this relation is linear. As we will see, this is the case in many practical situations such as the construction of corporate or sovereign bond yield curve, the construction of OIS discount curves, the construction of forward curves based on fixed-vs-floating interest-rate swaps or  the construction of implied default rates based on CDS spreads.\\

\subsection{Market-fit  and shape-preserving conditions}

%A term-structure curve is typically built by matching a set of market quotes $\mathbf{S}=(S_1, \ldots, S_n)$ corresponding to a series of financial products with increasing maturities $\mathbf{T}=(T_1, \ldots, T_n)$. The curve $t \rightarrow P(t_0, t)$ is compatible with the input set  $\mathbf{S}$ if for some payment time grid $(t_1, \ldots,t_{p_1}, \ldots, t_{p_n})$ with $t_{p_i} = T_i$, the column vector  $\mathbf{P} = (P(t_0, t_k))_{k=1, \ldots, p_n}$ is solution of the following rectangular linear system 
%\begin{equation}
%\label{eq:linear_sys}
%\mathbf{A}\cdot \mathbf{P}  = \mathbf{B}% \;\;\mbox{where}
%\end{equation}
%where
%%\begin{itemize}
%%\item $\mathbf{P}(\ba) = (P(t_0, t_k; \ba))_{k=1, \ldots, m}$ is the $m \times 1$ vector of elementary quantities that appear in the present value formula of instruments used to build the curve (\textcolor{blue}{see Examples 1 to 4}).\\
%$\mathbf{A}$ is a $n\times p_n$ matrix and $\mathbf{B}$ is a $n\times 1$ matrix. We moreover assume that, for any $i=1, \ldots, n$, the $i$-th line of the previous system corresponds to the market-fit condition of the financial product with maturity $T_i$. Note that matrices  $A$ and $B$ only depend on market quotes $\mathbf{S}$, on standard maturities  $\mathbf{T}$ and on products characteristics. For OIS discounting curve construction, $P = P^D$ and matrices $\mathbf{A}$ and  $\mathbf{B}$  can be easily extracted from system  (\ref{eq:OIS_market_fit}). For credit curve construction based on CDS spreads, $P = Q$ and matrices $\mathbf{A}$  and  $\mathbf{B}$ can be obtained from a discretized version of the system described by (\ref{eq:credit_system}).
%As a rectangular system ($n<p_n$),  (\ref{eq:linear_sys}) may admit several solutions\footnote{If maturity dates are strictly increasing, i.e., $T_1<\ldots<T_n$, $\mathbf{A}$ is a full rank matrix (rank $n$) and the solutions of  (\ref{eq:linear_sys}) evolves in a linear space with dimension equal to $n-p_n$.}.\\



The aim is to construct at some quotation date $t$ a term-structure function $T \rightarrow P(t, T)$, based on the observation  of a series of market quotes $S_1(t) \ldots, S_n(t)$ corresponding to the market value of $n$  financial instruments with time-to-maturities $T_1, \ldots, T_n$. In what follows, the quantity $T$ denotes a time length (as opposed to a calendar date), so that $P(t,T)$ corresponds to the value of the curve at time horizon $T$ or at calendar date $t+T$. The observation of market quotes at time $t$ provides a partial information on the curve at a set of time horizons or points $X =  \left(\tau_1, \ldots, \tau_m\right)$, i.e., at some calendar dates $t+\tau_1, \ldots, t+\tau_m$ which typically correspond to  payment dates of cash-flows.\\

The curve is (fully) compatible at time $t$ with market observations if the vector $P(t, X):=\left(P(t,\tau_1), \ldots, P(t,\tau_m)\right)^{\top}$ satisfies a linear system of the form 
\begin{equation}
\label{eq:market_fit}
A_t \cdot P(t, X)=\boldsymbol{b}_t,
\end{equation}
where $A_t$ is $n\times m$ real-valued matrix and $\boldsymbol{b}_t$ is a $n$-dimensional column vector. Of course, $A_t$ and $\boldsymbol{b}_t$ may depend on market quotes $S_1(t) \ldots, S_n(t)$, on the characteristics of the product cash-flows but also on the hypothesis made for assessing the value of these products at time $t$.
 Note that the number of observations $n$ is typically smaller than the number of points $m$, so that the solution of  system~\eqref{eq:market_fit} lives in a linear space  with dimension $m-n$.
\\

When constructing a financial term-structure, one may consider some additional information on the shape of the curve. For instance, the function $T \rightarrow P(t, T)$ may be known to be decreasing with respect to time horizon $T$ and its values may be bounded and belong to the interval $[0,1]$. This is  typically the case when one wants to construct a curve of discount factors (default-free zero-coupon bond prices) or an implied survival function (survival probabilities of a CDS reference entity). Violating these kinds of shape-preserving conditions results in term-structure functions that are typically not arbitrage-free.
%These shape-preserving conditions have to be considered in the process 
\\

\begin{rem}
\label{rem:surface}
The term-structure construction problem can be stated in a two-dimensional setting when the evolution of time (quotation dates) is added as a second dimension. In that case, the aim is to construct a surface $(t,T) \rightarrow P(t, T)$ based on a series of market quotes $S_1(t) \ldots, S_n(t)$ observed at several quotation dates $t=t_1, \ldots, t_N$. Note that the cash-flows characteristics of these products may depend on time $t$. In particular, the underlying maturity dates could be different at every time $t$. If for any date $t=t_1, \ldots, t_N$, the market-fit condition translates into a linear system, then this condition can again be expressed by a single system by concatenating  for every time $t=t_1, \ldots, t_N$ the $N$ systems  given as in Equation~\eqref{eq:market_fit}. The shape-preserving condition can be expressed as the intersection of shape-preserving conditions at time $t=t_1, \ldots, t_N$.
\end{rem}

% In what follows, $t_0$  denotes the quotation date. 
%of the products and we denote by $S$ the observed quotation for a product with a specific maturity time $T$.\\

%time, i.e., the date at which we observe the underlying product market quotes.\\ \\

In the following, we give some examples where  the construction of term-structures  involves combination of linear equality constraints and shape preserving conditions such as monotony or positivity.

\subsection{Classical examples of term-structures}

%\begin{ex}[Corporate or sovereign bond yield curve] 
%\label{ex:yield_curve}

In what follows, $t$  denotes a particular quotation date, i.e., the date at which the market quotes are observed for the underlying family of contingent products. 
For each example, we only consider a single financial product in this family and we provide the linear relation which characterize its value. These relation then corresponds to one particular row of the linear system~\eqref{eq:market_fit}.\\

\noindent \textbf{Corporate or sovereign bond yield curves}\\

Let $S$ be the observed market price  of a corporate or a sovereign bond with time-to-maturity  $T$ and with a fixed coupon rate $c$. The price $S$ and the coupon rate $c$ are expressed in percentage of invested nominal. The set of coupon payment dates is given by $(t+\ta_1,\ldots, t+\ta_p)$ where $\ta_1<\ldots<\ta_p=T$. The  year fraction $\delta_k$ represents the  time length $\ta_{k}- \ta_{k-1}$, $k=1, \ldots, p$ where  $\ta_0=0$. The present value of this bond can be defined as a linear combination of some default-free zero-coupon bonds, i.e., 
%$S$ can be decomposed as a linear combination of (fictitious 
\begin{equation}
\label{eq:bond}
 c\sum_{k=1}^{p}{\delta_k P^{B}(t,\ta_k)} + P^{B}(t, \tau_p) = S,
\end{equation}
where $P^{B}(t,\ta)$ represents the price at time $t$ of a default-free zero-coupon bond with time-to-maturity $\ta$. Note that, even if representation~\eqref{eq:bond} obviously relies on a default-free assumption, it is commonly employed as an intermediary step in the computation of the so-called bond yield-to-maturity.\footnote{The bond yield associated with time-to-maturity $T$ is defined as the constant rate of return $Y(t,T)$ such that the present value relation~\eqref{eq:bond} holds exactly when all the involved ZC bonds have this rate of return.} In this example, the curve $T\rightarrow P^B(t,T)$ will be inferred from a set of market fit conditions similar to~\eqref{eq:bond}, each of them corresponding to specific debt products but with different maturities.
As a result, this set of conditions can be easily represented in the form of a linear system as in Equation~\eqref{eq:market_fit}. In addition, the default-free assumption implies that the curve  $T\rightarrow P^B(t,T)$ is decreasing if arbitrage opportunities are precluded.\\
%\end{ex}


%\noindent \textbf{Example 1: Corporate or sovereign bond yield curve}\\
%
%Let $S$ be the observed market price (in percentage of nominal) of a corporate or a sovereign bond with maturity time $T$ and with a fixed coupon rate $c$. The set of coupon payment dates is given by $(t_1,\ldots, t_p)$ where $t_0<t_1<\ldots<t_p=T$ and $\delta_k$ represents the year fraction corresponding to period $(t_{k-1}, t_k)$, $
%k=1, \ldots, p$. The present value of this bond can be defined as a linear combination of some (fictitious) default-free zero-coupon bonds, i.e., 
%%$S$ can be decomposed as a linear combination of (fictitious 
%\begin{equation}
%\label{eq:bond}
% c\sum_{k=1}^{p}{\delta_k P^{B}(t_0,t_k)} + P^{B}(t_0, T) = S
%\end{equation}
%where $P^{B}(t_0,t)$ represents the price at time $t_0$ of a default-free zero-coupon bond with maturity $t$. Note that, even if the previous representation obviously relies on a default-free assumption, it is commonly used as an intermediary step in the computation of the so-called bond yield-to-maturity.\footnote{The bond yield for maturity time $T$ is defined as the constant rate of return $Y$ such that $P^{B}(t_0,t)$ has a return rate $Y$ for any time $t$, $t_0 \leq t \leq T$ and such that present value relation \ref{eq:bond} holds.}.\\
% 
% %and $\Delta_j$ the length of the period ($T_{j-1}, T_j)$ with $T_0=0$.\\ \\
%%, $c$ the coupon rate and $\Delta$ the cashflow frequency.  \\

\noindent \textbf{Discount curves based on Overnight Indexed Swaps par rates}\\

%\begin{ex}[Discounting curve based on OIS] 
%\label{ex:discount_curve}
Due to legal terms of standard collateral agreements, a possible choice to build discount curves is to use market quotes of OIS-like instruments \citep[see, for instance][for more details]{Hull2013}.
Let $S$ be the par swap rate of an overnight indexed swap with maturity $T$ and fixed leg payment dates $\ta_1<\ldots<\ta_p = T$. The  year fraction $\delta_k$ represents the  time length $\ta_{k}- \ta_{k-1}$, $k=1, \ldots, p$ where  $\ta_0=0$. For overnight-index swaps this time length is typically equal to one year.
The swap equilibrium relation takes the following linear form
\begin{equation}
\label{eq:OIS}
%S^{\mbox{OIS}} \sum_{j=1}^{n}  \Delta_j P^{\mbox{OIS}}(t,T_j)  + P^{\mbox{OIS}}(t, T_n) = 1
S \sum_{k=1}^{p}  \delta_k P^{D}(t,\ta_k)  = 1 - P^{D}(t, T), 
\end{equation}
where $P^{D}(t,\ta)$ is the discount factor at time $t$ associated with a time horizon $\tau$. In the previous equation, the left hand side represents the fixed leg present value whereas  the right hand side corresponds to the floating leg present value. For more details on the derivation of (\ref{eq:OIS}), the reader is referred to \citet{Fujii2010}.
\\

Under some circumstances, discount cuves can also be extracted from par rates of fixed versus Euribor swaps with Euribor tenor of 3 months or 6 months. This is the case for instance in the LTGA framework of Solvability 2  prudential regulation where the basis risk-free rates used for euro are constructed from the euro swap rates with a small adjustment for credit spread \citep[see, e.g.][]{CFO-CRO-2010}.
% \cite{van2014long})    
The resulting market-fit condition will exactly have the same form as~\eqref{eq:OIS}.\\

Note that, in the banking industry, discount curves are now understood as OIS based curves, see e.g. \cite{pallavicini2010interest}. The next example explains how to infer Euribor forward rates from quoted OIS rates and Euribor Swap rates.\\


%\end{ex}

%\noindent \textbf{Example 2: Discounting curve based on OIS}\\
%
%Due to the legal terms of standard collateral agreements, a possible choice to build discounting curves is to use quotations of OIS-like instruments.
%Let $S$ be the par swap rate of an overnight indexed swap with maturity $T$ and fixed leg payment dates $t_1<\ldots<t_p = T$. If for any $k=1, \ldots, p$, $\delta_k$ represents the year fraction corresponding to period $(t_{k-1}, t_k)$, 
%the swap equilibrium relation can be expressed in linear form with respect to some discount factors $P^{D}$, i.e., 
%
%\begin{equation}
%\label{eq:OIS}
%%S^{\mbox{OIS}} \sum_{j=1}^{n}  \Delta_j P^{\mbox{OIS}}(t,T_j)  + P^{\mbox{OIS}}(t, T_n) = 1
%S \sum_{k=1}^{p}  \delta_k P^{D}(t_0,t_k)  = 1 - P^{D}(t_0, T) 
%\end{equation}
%where $P^{D}(t_0,t)$ is the discount factor associated with maturity date $t$. In the previous equation, the left hand side represents the fixed leg present value whereas  the right hand side corresponds to the floating leg present value.
%\\


\noindent \textbf{Forward curves based on OIS and fixed versus Ibor-floating interest-rate swaps}\\

%\begin{ex}[Forward curve based on fixed-vs-Ibor-floating IRS] 

 Let $S$ be the observed par rate of an interest rate swap with maturity time $T$ and floating payments linked to a Libor or an Euribor rate associated with a tenor $j$ (typically, $j=3$ months or $j=6$ months). The fixed-leg payment scheme is given by $\ta_1<\cdots<\ta_p=T$ and  the floating-leg payment scheme is given by $\tilde{\ta}_1<\cdots<\tilde{\ta}_{q}=T$. 
  For most liquid products, payment on the fixed  leg are made with an annual frequency, so that $\ta_k$ corresponds to $k$ years ahead from the current date $t$. 
The  year fraction $\delta_k$ represents the  time length $\ta_{k}- \ta_{k-1}$, $k=1, \ldots, p$ ($\ta_0=0$) whereas the year fraction  $\tilde{\delta}_i$ represents the time length  $\tilde{\ta}_{i}-\tilde{\ta}_{i-1}$, $i=1, \ldots, q$ ($\tilde{\ta}_0=0$).   
Note that the length between two consecutive dates on the floating leg should correspond to the Libor or Euribor tenor, i.e., depending on the case, $\tilde{\delta}_i \simeq$ 3 months or $\tilde{\delta}_i \simeq$ 6 months. As a result, given an OIS discount curve $P^{D}$, the swap equilibrium relation can be represented in a linear form with respect to some forward Libor or Euribor rates, i.e.,  
\begin{equation}
\label{eq:IRS}
S \sum_{k=1}^{p}  \delta_k P^{D}(t,\ta_k) = \sum_{i=1}^{q} P^{D}(t, \tilde{\ta}_i) \tilde{\delta}_i F_j(t, \tilde{\ta}_{i} ),
\end{equation}
where $P^{D}(t, \tilde{\ta})$ is a risk-free discount factor at time $t$ for maturity $\ta$ and  $F_j(t, \tilde{\ta}_{i} ) := F(t,\tilde{\ta}_{i-1}, \tilde{\ta}_{i} )$ is the forward Libor or Euribor rate defined as the fixed rate to be exchanged at time $\tilde{\ta}_i$ against the $j$-tenor Libor or Euribor rate established at time $\tilde{\ta}_{i-1}$ so that the swap has zero value at time $t$. As in the previous example, the left hand side of~\ref{eq:IRS} represents the fixed leg present value whereas the right hand side corresponds to the floating leg present value. For more details, see, for instance \citet{Chibane2009}. Given a (pre-constructed) discount curve $P^D$ and  a set of swap par rates $S = S_1, \ldots, S_n$ corresponding to time-$t$ market quotes of Euribor or Libor swaps with maturities $T=T_1, \ldots, T_n$, the forward curve $T\rightarrow F_j(t, T)$ under construction has to satisfy a linear market-fit condition. This condition takes the form of a linear system whose each line is given by~\eqref{eq:IRS}. In addition, one can  require that forward rates are positive, so that the underlying pseudo zero-coupon prices form a decreasing function of time horizon. This particular shape preserving property (positivity) can also be enforced in the interpolation procedure we will described in the next section.\\


%\end{ex}


%\noindent \textbf{Example 3: Forward curve based on fixed-vs-Ibor-floating IRS}\\
%
% Let $S$ be the observed par rate of an interest rate swap with maturity time $T$ and floating payments linked to a Libor or an Euribor rate associated with a tenor $j$ (typically, $j=3$ months or $j=6$ months). The fixed-leg payment scheme is given by $t_1<\cdots<t_p=T$ and  the floating-leg payment scheme is given by $\tilde{t}_1<\cdots<\tilde{t}_{q}=T$. 
%  For most liquid products, payment on the fixed  leg are made with an annual frequency, so that $t_k$ is the business day corresponding to $k$ years after the current date $t_0$. The associated year fraction for the interval $(t_{k-1}, t_k)$ is denoted by $\delta_k$ whereas the year fraction for the interval $(\tilde{t}_{i-1}, \tilde{t}_{i})$ is denoted by $\tilde{\delta}_i$. Note that the length between two consecutive dates on the floating leg should be close to the Libor or Euribor tenor, i.e., $\tilde{\delta}_i \simeq j$. As a result, given a discounting curve $P^{D}$, the swap equilibrium relation can be represented in a linear form with respect to some forward Libor or Euribor rates, i.e.,  
%\begin{equation}
%\label{eq:IRS}
%S \sum_{k=1}^{p}  \delta_k P^{D}(t_0,t_k) = \sum_{i=1}^{q} P^{D}(t_0, \tilde{t}_i) \tilde{\delta}_i F_{j}(t_0, \tilde{t}_{i})
%\end{equation}
%where $P^{D}(t_0, t)$ is a risk-free discount factor at time $t_0$ for maturity $t$ and  $F_{\Delta}(t_0,t_{k})$ is the forward Libor or Euribor rate defined as the fixed rate to be exchanged at time $\tilde{t}_i$ against the $j$-tenor Libor or Euribor rate established at time $\tilde{t}_{i-1}$ so that the swap has zero value at time $t_0$. As in the previous example, the left hand side of \ref{eq:IRS} represents the fixed leg present value whereas the right hand side corresponds to the floating leg present value.\\ 
	

\noindent \textbf{Credit curves based on Credit Default Swaps spreads}\\
%\begin{ex}[Credit curve based on CDS] 
%\label{ex:credit_curve}

Let $S$ be the fair spread of a credit default swap with protection time  horizon $T$ and with premium payment dates $\ta_1<\cdots<\ta_p=T$. If we denote by $R$ the  expected recovery rate of the reference entity  and by 
$\delta_k$ the year fraction corresponding to time length $\ta_{k}- \ta_{k-1}$ ($\ta_0 = 0$), then 
%Discount factors will be denoted by $P$ the price of default free ZC prices whereas $\tilde{P}$ corresponds to (risk-neutral) survival probability of the underlying issuer. 
 the CDS swap equilibrium relation can be expressed as
 %in linear form with respect to some risk-neutral survival probabilities, i.e., 
% satisfies the following relation (under the assumption that the default payment would be executed at closest future premium payment date)
\begin{equation}
\label{eq:CDS}
  %S\sum_{j=1}^n \delta_j P(0,T_j) \tilde{P}(0,T_j) = (1-R)\sum_{j=1}^n P(0,T_j) \left(\tilde{P}(0,T_{j-1}) - \tilde{P}(0,T_{j})\right) 
  S \sum_{k=1}^p \delta_k P^D(t, \ta_k) Q(t, \ta_k)  = -(1-R)\int_{0}^T P^D(t, \ta) dQ(t,\ta),
\end{equation}
where $P^{D}(t, \ta)$ is the risk-free discount factor at time $t$ for time horizon $\ta$ and where $Q(t,\ta)$ is the probability (at time $t$) that the underlying reference entity has not defaulted before time horizon $\ta$. Then, $Q(t,\ta)$ is the survival probability of the  debt issuer in the time horizon $\ta$.  The left hand side of~\eqref{eq:CDS} represents the premium leg present value whereas the right hand side corresponds to the protection leg (or default leg) present value.
%probability at time-$t$ that the reference entity is still alive at time $T$. 
We implicitly assume here that recovery, default and interest rates are stochastically independent.
Using an integration by parts, it is straightfoward to show that survival probabilities $Q(t,\ta)$, $0 \leq \ta \leq T$,  are linked through the following linear relation~:
\begin{equation}
\label{eq:linear_CDS}
\begin{split}
S &\sum_{k=1}^p \delta_k P^D(t,\ta_k) Q(t, \ta_k) + (1-R) P^D(t,T)Q(t,T)\\ 
& \hspace{0.5 cm} + (1-R)\int_{0}^T f^D(t,\ta)P^D(t, \ta)Q(t,\ta)d\ta = 1 - R% \\
%& - (1-R)\left(1 - P^D(t,T)Q(t,T)\right)\\
%&  =  0
\end{split}
\end{equation}
%\begin{equation}
%\begin{aligned}
%& S \sum_{k=1}^n \delta_k P^D(t_0,t_k) Q(t_0, t_k) + (1-R) P^D(t_0,T)Q(t_0,T)\\ 
%& \hspace{0.5 cm} + (1-R)\int_{t_0}^T f^D(t_0,t)P^D(t_0, t)Q(t_0,t)dt = 1 - R% \\
%%& - (1-R)\left(1 - P^D(t,T)Q(t,T)\right)\\
%%&  =  0
%\end{aligned}
%\end{equation}
where $f^D(t,\ta)$ is the instantaneous forward rate\footnote{Instantaneous forward rates can be derived from discount factors through the following relation~: $f^D(t, \ta)P^D(t, \ta) = - \frac{\partial P}{\partial \ta} (t, \ta)$.} at time $t$ for  time horizon $\ta$.
In practice, the integral involved in the expression of the protection leg present value is classically discretized on the premium time grid $\ta_1<\cdots<\ta_n=T$, so that the continuous linear condition~\eqref{eq:linear_CDS} can be stated as a discrete one given by
\begin{equation}
\label{eq:discrete_linear_CDS}
\begin{split}
 &\sum_{k=1}^{p-1} \left(S\delta_k P^D(t,\ta_k) + (1-R)(P^D(t,\ta_{k-1}) - P^D(t,\ta_k))\right) Q(t, \ta_k) \\ 
& \hspace{0.5 cm} + \left(S\delta_p P^D(t,T) + (1-R) P^D(t,\ta_{p-1})\right) Q(t, T) = 1 - R.% \\
%& - (1-R)\left(1 - P^D(t,T)Q(t,T)\right)\\
%&  =  0
\end{split}
\end{equation}
At some fixed quotation date $t$, CDS protection is usually available for a set of liquid maturities $T_1, \ldots, T_n$. Then the construction of an implied survival function $T\rightarrow Q(t,T)$ consists in building a $[0,1]$-valued decreasing function that satisfies a system of $n$  linear equality constraints as of~\eqref{eq:discrete_linear_CDS}.


%\end{ex}


%\subsection{Literature}
%\subsection{Organization of the paper}

%\subsection{Matrix formulation of the market-fit condition}
%
%The objectif is to construct a term-structure function $(t,T) \rightarrow P(t, T)$ based on the observation at any time $t=t_1, \ldots, t_N$ of  a series of market quotes $S_1(t) \ldots, S_n(t)$ corresponding to the market value of $n$ benchmark financial instruments. We assume that, for any quotation dates $t=t_1, \ldots, t_N$, this set of products generate a stream of cash-flows at dates $\tau_j$, $j=1, \ldots, m$. 
%\begin{equation*}
%A_t \cdot P(t, X)=\boldsymbol{b}_t, \qquad, \ t=t^{(1)}, \ldots,t^{(T)}
%\end{equation*}
%where $P(t, X) = \left(P(t, \tau_1), \ldots, P(t,\tau_m)\right)$\\
%
%
%
%Note that, term-structure surface can also be constructed based on the information on market quotes collected at different time $t=t_1, \ldots, t_N$. In that case, one can translate the series of linear constraints into a single system of the form : concatenation of the previous $n$ systems...


%\section{Motivation of the constraints}
%\label{Sec:motiv_constraints}
%
%Here are the references mentioned by the reviewers.
%\begin{itemize}
%\item Reviewer 2 : \cite{LauriniOhashi2015} (cité dans Nugget), \cite{Fengler:2015} (cité), \cite{Laurini:2010} (cité)
%\item Reviewer 3 : \cite{Kanevski2008} (cité à la fin), \cite{Benth:2015} (cité à la fin), \cite{Steeley:2008} (cité à la fin), \cite{nelson1987parsimonious}, \cite{Svensson:1994}
%\end{itemize}

\section{Kriging under linear equality constraints}

\label{IntModels}

The term-structure construction approach we propose relies on an extension of kriging to linear equality and shape preserving constraints. In this section, we give a formal presentation of this interpolation technique when only linear equality constraints are considered. Section~\ref{IntModIConst} explains how this technique can be generalized when some monotonicity constraints are added.\\ 

Kriging or Gaussian process regression is a method of interpolation for which the interpolated values are modeled by a Gaussian Process (GP) with a prior covariance function. This method is widely used in the domain of spatial analysis and computer experiments \citep[see, e.g.][]{Rasmussen:2005:GPM:1162254}. More formally, we consider the model $y=f(\boldsymbol{x})$, where $f$ is an unknown real-valued function of a $d$-dimensional input variable $\boldsymbol{x}\in\mathbb{R}^d$. In the case where the computation of $f$ is expensive and time-consuming or in the case where $f$ is known only at some input locations, this function can be estimated using so-called \emph{kriging} techniques. 
In that case, $f$  is seen as a realization of a Gaussian Process (GP) $Y$ defined as
\begin{equation*}
Y(\boldsymbol{x}):=\mu(\boldsymbol{x})+Z(\boldsymbol{x}),
\end{equation*}
where the deterministic function $\mu \ : \ \boldsymbol{x}\in \mathbb{R}^d \longrightarrow \mu(\boldsymbol{x})\in\mathbb{R}$
is the mean of $Y$ and $Z$ is a zero-mean GP with covariance function
\begin{equation*}
K \ : \ (\boldsymbol{x},\boldsymbol{x}')\in \mathbb{R}^d\times\mathbb{R}^d \longrightarrow K(\boldsymbol{x},\boldsymbol{x}')=\cov(Y(\boldsymbol{x}),Y(\boldsymbol{x}'))\in\mathbb{R}.
\end{equation*}
The quantity $K(\boldsymbol{x},\boldsymbol{x}')$ is thus the covariance of the values of the random field $Y$ at two input locations $\boldsymbol{x}$ and $\boldsymbol{x}'$.
We assume that the covariance function $K$ is such that the random field $Y$ have continuous and differentiable sample paths with probability one, see \cite{abrahamsen1997review}. Figure~\ref{GPsimkernels} (right) gives an illustration of some  Gaussian processes sample paths.  In numerical illustrations of Section~\ref{EmpInv}, we consider Gaussian processes with $d$-dimensional covariance functions given as a tensor product, i.e., for $\boldsymbol{x}=(x_1,\ldots,x_d)$ and $\boldsymbol{x}'=(x'_1,\ldots,x'_d)$~:
\begin{equation*}
K(\boldsymbol{x},\boldsymbol{x}')=\sigma^2\prod_{i=1}^dC_i(x_i-x'_i,\theta_i),
\end{equation*}
where $\boldsymbol{\theta}=(\theta_1,\ldots,\theta_d)\in\mathbb{R}^d$ and $\sigma^2$ are respectively called \textit{length} and \textit{variance hyper parameters}.
The function $C_i$ are kernel correlation functions which depend on the length parameter $\theta_i$ and on $x_i-x'_i, \ i=1,\ldots,d$, see Table~\ref{kernel} for some popular kernel correlation functions. Note that the length parameter $\theta$ can be interpreted as a correlation parameter\footnote{For a Gaussian covariance kernel, it can be shown that, for any $n$,  increasing $\boldsymbol{\theta}$ yields an increase of the vector $(Y(x^{1}), \ldots, Y(x^{n}))$ with respect to the supermodular order, a stochastic order which is well-known for comparing the degree of dependence.} as it controls the degree of dependence amongst the values of the Gaussian process at any two points. The parameter $\sigma$ controls the initial Gaussian process variance.
% \hassan{in one dimension}.

\subsection{Classical kriging}
In classical kriging, the real function $f$ is known to take some values $y_1,\ldots,y_n$ at some $d$-dimensional design points $x^{(1)},\ldots,x^{(n)}$, so that $f(X)=\boldsymbol{y}$, where the design points are given as the rows of the $n\times d$ matrix $X=\left(x^{(1)},\ldots,x^{(n)}\right)^\top\in\mathbb{R}^{n\times d}$, $f(X)=\left(f(x^{(1)}),\ldots,f(x^{(n)})\right)^\top\in\mathbb{R}^n$ and $\boldsymbol{y}=(y_1,\ldots,y_n)^\top\in\mathbb{R}^n$. One advantage of using a GP emulator is that, conditionally to the observation data $\boldsymbol{y}$, the conditional process $Y \suchthat Y(X)=\boldsymbol{y}$ is still a GP. This process is characterized by its (marginal) mean 
\begin{equation}
\label{eq:kriging_mean}
\eta(\boldsymbol{x})=\mu(\boldsymbol{x})+\boldsymbol{k}(\boldsymbol{x})^\top\mathbb{K}^{-1}(\boldsymbol{y}-\boldsymbol{\mu}), \quad  
\boldsymbol{x}\in\mathbb{R}^d
\end{equation}
and its covariance function $\tilde{K}$ given by
\begin{equation}
\tilde{K}(\boldsymbol{x},\boldsymbol{x}')= K(\boldsymbol{x},\boldsymbol{x}')-\boldsymbol{k}(\boldsymbol{x})^\top\mathbb{K}^{-1}\boldsymbol{k}(\boldsymbol{x}'),\quad  
\boldsymbol{x},\boldsymbol{x}' \in\mathbb{R}^d
\end{equation}
%marginal distributions of $Y(\boldsymbol{x})$ are given by, for all $\boldsymbol{x}\in\mathbb{R}^d$~:
%\begin{equation*}
%Y(\boldsymbol{x})\mid Y(X)=\boldsymbol{y} \ \sim \ \mathcal{N}(\eta(\boldsymbol{x}),\tau^2(\boldsymbol{x})), 
%\end{equation*} 
%where
%\begin{eqnarray*}\left\{\begin{array}{ll}
%\eta(\boldsymbol{x})=\mu(\boldsymbol{x})+\boldsymbol{k}(\boldsymbol{x})^\top\mathbb{K}^{-1}(\boldsymbol{y}-\boldsymbol{\mu})\\
%\tau^2(\boldsymbol{x})=K(\boldsymbol{x},\boldsymbol{x})-\boldsymbol{k}(\boldsymbol{x})^\top\mathbb{K}^{-1}\boldsymbol{k}(\boldsymbol{x})
%\end{array}\right.
%\end{eqnarray*}
where $\boldsymbol{\mu}=\mu(X)=\left(\mu(x^{(1)}),\ldots,\mu(x^{(n)})\right)^\top\in\mathbb{R}^n$ is the trend vector at the design points, $\mathbb{K}$ is the covariance matrix
of $Y(X)$ and $\boldsymbol{k}(\boldsymbol{x})=\left(K\left(\boldsymbol{x},x^{(1)}\right),\ldots,K\left(\boldsymbol{x},x^{(n)}\right)\right)^\top$ is the vector of covariance between $Y(\boldsymbol{x})$ and $Y(X)$. The conditional mean $\eta(\boldsymbol{x})$ given the observation data $Y(X)=\boldsymbol{y}$ is the Best Linear Unbiased Estimator (BLUE) of $Y(\boldsymbol{x})$, which is known as \textit{kriging mean} \citep[see][]{jones1998}. 
%The conditional variance $\tau^2(\boldsymbol{x})$ is called the \textit{kriging variance}. 
One remarkable property is that the covariance function $\tilde{K}$ of the conditional Gaussian process does not depend on the observation data $\boldsymbol{y}$. 
%\hassan{Notice that for characterizing the conditional Gaussian process, one also need the conditional covariance function  
%\begin{equation*}
%\cov(Y(\boldsymbol{x}),Y(\boldsymbol{x}')\suchthat Y(X)=\boldsymbol{y})=K(\boldsymbol{x},\boldsymbol{x}')-\boldsymbol{k}(\boldsymbol{x})^\top\mathbb{K}^{-1}\boldsymbol{k}(\boldsymbol{x}').
%\end{equation*} }
In addition, the regularity of the kriging mean predictor function $\eta$ inherits from  the regularity of the mean function $\mu$ and from the regularity of the
covariance function $K$ of the original GP $Y$. %which is symmetric and positive definite. 
Then, the choice of this covariance function  is essential because it drives the smoothness of the kriging metamodel. Table~\ref{kernel} gives some popular  kernel correlation functions, ordered by decreasing degree of smoothness.
 %their corresponding degree of regularity. 
 %like Gaussian function, Mat\'ern function (with parameters $\lambda=3/2,\lambda=5/2,\cdots$) and exponential function (which is a Mat\'ern function with parameter $\lambda=1/2$). 
%Notice that these functions are placed in order of smoothness, the Gaussian function is smooth function and the exponential function is of class \hassan{$\mathcal{C}^0$ (continuous but non-derivable functions)} \cite{Rasmussen:2005:GPM:1162254}. 
Figure~\ref{GPsimkernels} shows three alternative covariance functions with  associated Gaussian process sample paths.  
\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GPkernels.pdf}
%\subfloat[]{\label{GPkernels}\includegraphics[scale=.4]{figures/GPkernels.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GPsimulations.pdf}
%\subfloat[]{\label{GPsimulations}\includegraphics[scale=.4]{figures/GPsimulations.pdf}}
\end{minipage}
\caption{Some Gaussian process covariance functions (left) and associated sample paths (right).
The covariance parameters are fixed to $(\theta, \sigma)=(0.3,\, 1.0)$.}
\label{GPsimkernels}
\end{figure}

  
\begin{table}[hptb]
\centering
\caption{Some popular kernel correlation functions $C(x-x',\theta)$ used in kriging methods.}
\label{kernel}
\begin{tabular}{ccc}
\hline \hline
Name  &  Expression &  Class   \\   
\hline
Gaussian       & $\exp\left(-\frac{(x-x')^2}{2\theta^2}\right)$ 	& $\mathcal{C}^{\infty}$       \\
Mat\'ern 5/2   & $\left(1+\frac{\sqrt{5}| x-x' |}{\theta}+\frac{5(x-x')^2}{3\theta^2}\right)\exp\left(-\frac{\sqrt{5} | x-x'|}{\theta}\right)$
& $\mathcal{C}^2$      \\
Mat\'ern 3/2   &  $\left(1+\frac{\sqrt{3}| x-x' |}{\theta}\right)\exp\left(-\frac{\sqrt{3}|x-x'|}{\theta}\right)$   & $\mathcal{C}^1$      \\
Exponential    &  $\exp\left(-\frac{|x-x'|}{\theta}\right)$ & $\mathcal{C}^0$      \\   
\hline
\end{tabular}
\end{table}


\subsection{Extension to linear equality constraints}\label{ILC}

The previous setting can be generalized by considering linear equality constraints instead of pure interpolation constraints. This is of primary importance if one wants to construct term-structures which are compatible with linear market-fit conditions. We are under the situation where $n$ relevant financial products are considered to construct the curve and their market quotes provide information on the curve at the $m$ points $x^{(1)},\ldots,x^{(m)}$. Then, the (unknown) real function $f$ satisfies some linear constraints of the form
\begin{equation}\label{initLinConst}
A\cdot f(X)=\boldsymbol{b},
\end{equation}
where $A$ is a given matrix of dimension $n\times m$, $n,m\in\mathbb{N}$, $X=\left(x^{(1)},\ldots,x^{(m)}\right)^\top\in\mathbb{R}^{m\times d}$, $f(X)=\left(f(x^{(1)}),\ldots,f(x^{(m)})\right)^\top \in \mathbb{R}^m$ and $\boldsymbol{b}\in\mathbb{R}^n$. 
%\ares{In classical kriging, the real function $f$ is known to take some values $y_1,\ldots,y_n$ at some $d$-dimensional design points $x^{(1)},\ldots,x^{(n)}$, so that $f(X)=\boldsymbol{y}$, where the design points are given as the rows of the $n\times d$ matrix $X=\left(x^{(1)},\ldots,x^{(n)}\right)^\top\in\mathbb{R}^{n\times d}$, $f(X)=\left(f(x^{(1)}),\ldots,f(x^{(n)})\right)^\top\in\mathbb{R}^n$ and $\boldsymbol{y}=(y_1,\ldots,y_n)^\top\in\mathbb{R}^n$}
%Note that in practice, $m$ is greater than $n$ so that the linear system admits 
%In the following, we focus on this type of linear equality constraints which generalizes the above equality ones. 
In that case, the conditional process $Y \suchthat A\cdot Y(X)=\boldsymbol{b}$ is still a Gaussian process with mean
\begin{equation}
\label{eq:mean}
\eta(\boldsymbol{x})=\mu(\boldsymbol{x})+\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}(\boldsymbol{b}-A\boldsymbol{\mu}), \quad  
\boldsymbol{x}\in\mathbb{R}^d
\end{equation}
and covariance function %$\tilde{K}$ such that 
\begin{equation}
\label{eq:cov}
\tilde{K}(\boldsymbol{x},\boldsymbol{x}')= K(\boldsymbol{x},\boldsymbol{x}')-\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}A\boldsymbol{k}(\boldsymbol{x}'),\quad  
\boldsymbol{x},\boldsymbol{x}' \in\mathbb{R}^d
\end{equation}
where $\boldsymbol{\mu}=\mu(X)=\left(\mu(x^{(1)}),\ldots,\mu(x^{(m)})\right)^\top\in\mathbb{R}^m$ is the trend vector at the design points, $\mathbb{K}$ is the covariance matrix
of $Y(X)$ and $\boldsymbol{k}(\boldsymbol{x})=\left(K\left(\boldsymbol{x},x^{(1)}\right),\ldots,K\left(\boldsymbol{x},x^{(m)}\right)\right)^\top$ is the vector of covariance between $Y(\boldsymbol{x})$ and $Y(X)$.
%$$
%\did{\cov(Y(\boldsymbol{x}), Y(\boldsymbol{x}') \mid A\cdot Y(X)=\boldsymbol{b})}= K(\boldsymbol{x},\boldsymbol{x}')-\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}A\boldsymbol{k}(\boldsymbol{x}'),\quad  
%\boldsymbol{x},\boldsymbol{x}' \in\mathbb{R}^d
%$$
%$$
%\cov(\boldsymbol{x}, \boldsymbol{x}')=K(\boldsymbol{x},\boldsymbol{x}')-\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}A\boldsymbol{k}(\boldsymbol{x}'),\quad  
%\boldsymbol{x},\boldsymbol{x}' \in\mathbb{R}^d
%$$
%the marginal distribution of the conditional GP is given as, for any $\boldsymbol{x}\in\mathbb{R}^d$  
%\begin{equation*}
%Y(\boldsymbol{x})\mid A\cdot Y(X)=\boldsymbol{b} \ \sim \ \mathcal{N}(\eta(\boldsymbol{x}),\tau^2(\boldsymbol{x})), 
%\end{equation*} 
%where
%\begin{eqnarray*}\left\{\begin{array}{ll}
%\eta(\boldsymbol{x})=\mu(\boldsymbol{x})+\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}(\boldsymbol{b}-\boldsymbol{\mu})\\
%\tau^2(\boldsymbol{x})=K(\boldsymbol{x},\boldsymbol{x})-\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top\right)^{-1}A\boldsymbol{k}(\boldsymbol{x})
%\end{array}\right.
%\end{eqnarray*}
Note that when $A$ is a square identity matrix, the linear constraints become interpolation constraints.

\subsection{Kriging in the presence of  noisy observations}\label{subsec:noisy}

The previous framework can be further extended to situations where market observations cannot be considered to be fully reliable due to, e.g., market microstructure effect. We assume that this uncertainty blurred the market information on $f$ in such a way that the previous measurement equation \eqref{initLinConst} is flawed by an additive error term $\epsilon$ :
$$
 \boldsymbol{b} = A\cdot f(X) + \epsilon.
$$
We consider that this error term is one realization of an independent zero-mean Gaussian noise $\boldsymbol{\varepsilon}$ in $\R^{m}$ with  covariance matrix $\Sigma$. 
%For instance, when  $\boldsymbol{b}$ represents mid-price of some reference products in a market situation where bid-ask spreads are large, $\Sigma$ can be chosen in such a way that  standard deviation of   $\boldsymbol{\varepsilon}$ is close to the differences between ask and mid prices.
%This is the case when market observations are not representative of the reference product fundamental values due to market microstructure contamination. For instance, when the bid-ask spreads are large, using mid-price to construct term-structure  function may result in flawed market fit condition. We then consider that we only have access to the transformation
%$A\cdot f(X)$ of $f(X)$ through an approximation  $\boldsymbol{b}  = A\cdot f(X) +  \boldsymbol{\epsilon}$ of it. When it is reasonable to assume that $\boldsymbol{\epsilon}$ is one realization of a “noise” random variable $\boldsymbol{\varepsilon}$, it is still possible to derive kriging approximations of function $f$. %For instance, one can consider that the noise $\varepsilon$ has a 
%Let us first assume that  market quotes appear in the right hand side of the linear equality constraint as it is the case for \eqref{eq:bond}. 
%This issue can be tackled by assuming that the value of the  function $f$ at point $X = ()$ satisfy a linear system of the form
%$$
% \boldsymbol{b}  = A\cdot f(X) +  \boldsymbol{\varepsilon}
%$$
%where $\boldsymbol{\varepsilon}$ is an independent zero-mean gaussian vector in $\R^{m}$ and with  covariance matrix $\Sigma$.   
%When  $\boldsymbol{b}$ is a vector of market quotes as it is the case in \eqref{eq:bond}, $\Sigma$ can be chosen as a diagonal matrix with diagonal  term $(\sigma^{(1)}, \ldots, \sigma^{(m)})$ such that $\sqrt{\sigma^{(j)}}$ is close to the difference between ask price and mid price.\\ 
%One can consider that observed quotes only represent an approximation of the true intrinsic value of the reference products. 
%We can then consider that, for a given reference product, a noisy realization of  market quote $S+\epsilon$ is observed instead of the `true' value $S$. 
%When it is reasonable to assume that $\epsilon$ is a realization of a “noise” random variable $\varepsilon$, it is still possible to derive financial term-structures as kriging approximations. Let us first assume that  market quotes appear in the right hand side of the linear equality constraint as it is the case for \eqref{eq:bond}. In this situation, 
%Each term of the diagonale corresponds to the uncertainty attached to the market fit condition. 
If the measurement error $\boldsymbol{\varepsilon}$ is assumed to be  independent of the Gaussian process $Y$,  the conditional process $Y \suchthat  \boldsymbol{b} = A\cdot Y(X) + \boldsymbol{\varepsilon}$ is still a Gaussian process with mean
\begin{equation}
\label{eq:mean_noisy}
\eta(\boldsymbol{x})=\mu(\boldsymbol{x})+\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top + \Sigma\right)^{-1}(\boldsymbol{b}-A\boldsymbol{\mu}), \quad  
\boldsymbol{x}\in\mathbb{R}^d
\end{equation}
and covariance function %$\tilde{K}$ such that
\begin{equation}
\label{eq:cov_noisy}
\tilde{K}(\boldsymbol{x},\boldsymbol{x}')= K(\boldsymbol{x},\boldsymbol{x}')-\left(A\boldsymbol{k}(\boldsymbol{x})\right)^\top \left(A\mathbb{K}A^\top + \Sigma\right)^{-1}A\boldsymbol{k}(\boldsymbol{x}'),\quad  
\boldsymbol{x},\boldsymbol{x}' \in\mathbb{R}^d
\end{equation}
where $\boldsymbol{\mu}=\mu(X)=\left(\mu(x^{(1)}),\ldots,\mu(x^{(m)})\right)^\top\in\mathbb{R}^m$ is the trend vector at the design points, $\mathbb{K}$ is the covariance matrix
of $Y(X)$ and $\boldsymbol{k}(\boldsymbol{x})=\left(K\left(\boldsymbol{x},x^{(1)}\right),\ldots,K\left(\boldsymbol{x},x^{(m)}\right)\right)^\top$ is the vector of covariance between $Y(\boldsymbol{x})$ and $Y(X)$.
Note that the only difference compared to noiseless kriging equations \eqref{eq:mean} and  \eqref{eq:cov} is the replacement, at every occurrence, of the covariance matrix $A\mathbb{K}A^\top$ by an inflated matrix $A\mathbb{K}A^\top + \Sigma$. In addition, the curve associated with the kriging mean function \eqref{eq:mean_noisy} does not satisfy the noiseless market fit condition.\\



%e.g. Márcio Poletti Laurini, Alberto Ohashi, A noisy principal component analysis for forward rate curves, European Journal of Operational Research, Volume 246, Issue 1, 1 October 2015, Pages 140-153.

The possibility to handle the presence of noisy observations is quite important. The presence of measurement errors and observations with contaminations may have important consequences on the dynamics of yield curves, as noticed for example in~\cite{LauriniOhashi2015}. The previous framework can be used to construct term structure functions in the presence of illiquid securities. Assume that the vector $\boldsymbol{b}$ represents mid-prices of some reference products in a market situation where bid-ask spreads are large. For instance, the covariance matrix $\Sigma$ can be chosen in such a way that the range of $(\boldsymbol{b}-\boldsymbol{\varepsilon})$ is concentrated on the bid-ask interval.
One possible way to do that is to define $\Sigma$ as a diagonal matrix where each positive term is the square difference between ask and mid prices.\\ 

Note moreover that this extension can also be useful when one wants to explicitly look for an attractive trade-off between precision and curve smoothness. Indeed, considering that observed information is uncertain allows to weaken the market-fit condition which in turn benefits to curve smoothness.
%This is the case in situations where market quotes do not need to be matched perfectly.
Our methodology can be then adapted both for best-fitting of noisy market prices and for the construction of an exact interpolatory term structure to a set of liquid instruments.

 %\ares{To be detailed}.

%REVIEWER 2: 1 - In page 3 the authors write "In addition, kriging can also handle the presence of noisy observations (using the so-called nugget effect)". This feature of the methodology is quite important, since the presence of measurement errors and observations with contaminations may have important consequences on the dynamics of yield curves, e.g. Márcio Poletti Laurini, Alberto Ohashi, A noisy principal component analysis for forward rate curves, European Journal of Operational Research, Volume 246, Issue 1, 1 October 2015, Pages 140-153. A simulation study showing the properties of the method in the presence of yield curve contaminations would be an important addition to the work.






% subject to uncertainty due, for instance, to the lack of liquidity of the underlying security market.  correct values $S_{1}, \ldots, S_{m}$ of the reference products at maturity dates $x^{(1)}, \ldots, x^{(m)}$ are observed with a 

%Laurini et Ohashi 2014
%It can also be induced by observational errors due to market microstructure effects. The existence of an underlying bid and ask bond price structure contaminates the yield and forward rate curves (see e.g Mizrach and Neely (2011) and Goyenko et al. (2011)).

% This is the case when one wants to include illiquid securities whose price does not need to be matched perfectly or one wants to explicitly look for an attractive trade-off between precision and curve smoothness.
%
% Indeed, when the liquidity of the underlying reference products is not sufficient, could we incorporate bid-ask spreads in the curve construction method ?
%
% Following the convenient particular case of Monte Carlo simulations, we finally make the assumption that εi ∼ N(0, τi2) (1 ≤ i ≤ n) independently. Note that although not explicitly addressed in this paper, the case of multi-Gaussian vector of εi’s with prescribed covariance matrix is a straightforward generalization of the model




\section{Kriging under additional monotonicity constraints}\label{IntModIConst}

As mentioned in Section~\ref{Sec:term-structure}, the studied real function $f$ may be known to satisfy some shape-preserving constraints such as monotonicity or positivity. 
\emph{No-Arbitrage constraints} lead to a theoretical need of monotonicity in curve constructions. For example, quantities such as prices of default-free zero-coupon bond, discount factors or implied survival probabilities shall be non-increasing with respect to time horizons under the no-arbitrage condition.  In this section, we propose to extend constrained spline techniques to constrained kriging in order to build term-structure functions.\\ 

The use of monotonic splines  for term-structure construction has been a subject of great interest in the literature. This has been investigated by among other \cite{Barzanti-Corradi99}, \cite{ramponi2003adaptive},  \cite{Chiu_el2008}. The relationship of shape and monotonicity restrictions with no-arbitrage conditions is detailed for example in the recent literature, see~\cite{Laurini:2010} and~\cite{Fengler:2015}. 
%Our setting thus corresponds to \emph{classical constrained setting}, which appears e.g. in the constrained splines literature.
Additional motivations of the use of monotonicity constraints will be developed in the numerical Section~\ref{sec:motivMonotonie}.\\

%However, classical kriging techniques do not consider monotonicity constraints. 
A natural extension of kriging under monotonicity constraints is to consider conditional Gaussian processes with monotone paths. However, the difficulty is that the conditional monotone process is not a Gaussian process any more. We adopt here the approach introduced by \cite{maatouk:hal-01096751} where Gaussian processes are approximated by finite-dimensional versions, so that the monotonicity constraints can  be checked very efficiently in the entire domain. \\ 

In the following, we consider linear equality constraints as described in Section~\ref{ILC} and we explain how to incorporate supplementary monotonicity constraints. Note that, even if we focus on monotonicity constraints, other shape-preserving constraints can be incorporated using similar ideas \citep[see][]{maatouk:hal-01096751}.
%\did{We motivate the choice of using specific monotone interpolation techniques instead classical non constraint ones in Section~\ref{sec:motivMonotonie}. We first introduce monotonicity constraints under }
We first introduce monotonicity constraints under 
 the  one-dimensional case, where one has to retrieve a monotonic curve at a given quotation date (Section~\ref{ODC}). We then explain how to extend the one-dimensional kriging construction method to dimension two, where information at several quotation dates can be jointly used  (Section~\ref{TDC}).% to construct a surface of term-structures (Section~\ref{TDC}).
%As mentioned above, the real function (physical phenomena) may be known to satisfy some inequality constraints
%such as monotonicity. For example, discount factors take values between $0$ and $1$ and they are known to be non-increasing with respect to maturities. The difficulty is that the usual kriging model does not respect monotonicity constraints. The solution is to incorporate monotonicity constraints into a GP emulator, which is also a hard problem as the resulting conditional process is not a GP. In this framework, a new model to incorporate monotonicity constraints into a GP emulator is presented in \cite{maatouk:hal-01096751}. 

%The inequality constraints are forced by constructions and respected in the whole input domain. 





\subsection{One dimensional case}\label{ODC}

In this section, we assume that the input variable $x$ belongs to an interval $\Dom = [\underline{x}, \overline{x}]$  of $\mathbb{R}$ and we consider an original Gaussian process $Y$ with covariance function $K$. For simplicity, we assume that $Y$ is a zero-mean GP. 
The aim of this section is to explain how to construct a process that both satisfies linear equality constraints  and monotonicity constraints.\\

 The first step is to approximate the original Gaussian process $Y$ by a finite-dimensional version which is monotonic under some simple finite-dimensional linear inequality conditions.
  We begin by discretizing the input interval $\Dom$ into  a regular subdivision $u_0<\ldots < u_N$ with $u_0=\underline{x}$, $u_N=\overline{x}$ and with a constant mesh $\delta$, so that $u_j=u_0+j\delta$, $j=0, \ldots, N$. We then consider an associated set of basis functions  $\phi_j$,  $j=0, \ldots, N$ defined as 
\begin{equation*}
\phi_j(x)=\int_{\underline{x}}^xh_j(u)du, \qquad x\in \Dom,
\end{equation*}
where $h_j(x):=\max\left(1-\frac{|x-u_j|}{\delta},0\right)$ is a hat function centered at the $j^{\text{th}}$ knot $u_j$ of the input subdivision.  Note that the basis functions  $\phi_j$,  $j=0, \ldots, N$ are increasing as primitives of positive functions.
%where the basis functions are chosen as primitive of hat functions~: 
We then define the finite-dimensional approximation of $Y$ on $\Dom$ as the process $Y^N$ such that 
%  
%  (see \cite{maatouk:hal-01096751})
% 
% The main idea in \cite{maatouk:hal-01096751} is to approximate the original zero-mean GP $Y$ by a finite-dimensional one of the form~:
\begin{equation}\label{propModel}
Y^N(x)=\eta+\sum_{j=0}^N\xi_j\phi_j(x),\qquad x\in\Dom,
\end{equation} 
where $\boldsymbol{\xi}=(\eta,\xi_0,\ldots,\xi_N)^\top$ is a zero-mean Gaussian vector.
If the Gaussian process $Y$ has almost surely  differentiable paths, choosing $\eta = Y(u_0)$ and $\xi_{j} = Y'(u_j)$, $j=0,\ldots,N$ guarantees that the finite-dimensional process $Y^N$ uniformly converges on $D$ towards $Y$ almost surely as $N$ tends to infinity \citep[see][]{maatouk:hal-01096751}. In that case, the
% the uniform convergence of $Y^N$ towards $Y$ is stated in \cite{maatouk:hal-01096751}. 
 covariance matrix $\Gamma^N$  of $\boldsymbol{\xi}$ is given as
\begin{equation}
\label{eq:cov_mat_one_dim}
\Gamma^N=\left[\begin{matrix}
K(u_0,u_0) && \frac{\partial K}{\partial x'}(u_0,u_j) \\\\
\frac{\partial K}{\partial x}(u_i,u_0) && \frac{\partial^2K}{\partial x\partial x'}(u_i,u_j)\\
\end{matrix}
\right]_{0\leq i,j\leq N},
\end{equation}
where $K$ is the covariance matrix of the original GP $Y$ and $u_j, \ j=0,\ldots,N$ are the knots of the input subdivision. 
%When $\xi_{j} = Y(u_j)$, $j=0,\ldots,N$,  the uniform convergence of $Y^N$ towards $Y$ is stated in \cite{maatouk:hal-01096751}. It can also be shown that $Y^N$ is again a Gaussian process and that its covariance function $K_N$ uniformly converges to the covariance function $K$ of $Y$. As can be seen in the next proposition, the monotonicity constraint reduces to a finite-dimensional (linear) inequality constraints on the underlying Gaussian vector $\boldsymbol{\xi}$.

\begin{pro}[Monotonicity]
\label{pro:monotone_1D}
The process $Y^N$  defined in Equation~\eqref{propModel} is non-increasing (resp. non-decreasing) \textit{if and only if} all the coefficients $\xi_j, \ j=0,\ldots,N$ are nonpositive (resp. nonnegative).
\end{pro}

\begin{proof}
If $\xi_j, \ j=0,\ldots,N$ are nonpositive then, since $\phi_j$ are increasing, then $Y^N$ is increasing. For the reverse implication, let us first notice that
the derivative of the basis functions $\phi_j, \ j=0,\ldots,N$ are such that
\begin{eqnarray*}
\phi_j'(u_k)=h_j(u_k)=\indic{j=k}=\left\{\begin{array}{ll}
1 & \mbox{if} \ j=k\\
0 & \mbox{if} \ j\neq k 
\end{array}\right.
\end{eqnarray*}
Thus, the derivative of the process $Y^N$ at any knots $u_k, \ k=0,\ldots,N$ is
\begin{equation*}
\left(Y^N\right)'(u_k)=\sum_{j=0}^N\xi_j\phi'_j(u_k)=\xi_k,
\end{equation*}
which concludes the proof.
%Then, the model defined in (\ref{propModel}) can be reformulated as
%\begin{equation*}
%Y^N(x)=Y^N(0)+\sum_{j=0}^N(Y^N)'(u_j)\phi_j(x),
%\end{equation*}
%since $\phi_j(0)=0, \ j=0,\ldots,N$. 
\end{proof}
The choice of the basis functions $\phi_j$ and of $\Gamma^N$ depends on the type of shape-preserving constraints. Other type of basis function can be used for other constraints  \citep[for more details, see][]{maatouk:hal-01096751}.\\


%\paragraph{Convergence.} Let us recall that the basis functions $\phi_j, \ j=0,\ldots,N$ verify~: 
%\begin{eqnarray*}
%\phi_j'(u_k)=h_j(u_k)=\indic{j=k}=\left\{\begin{array}{ll}
%1 & \mbox{if} \ j=k\\
%0 & \mbox{if} \ j\neq k 
%\end{array}\right.
%\end{eqnarray*}
%where $\phi'$ is the derivative of the basis function. Now, if we take the derivative of the Gaussian process $Y^N$ at any knots $u_k, \ k=0,\ldots,N$, we obtain 
%\begin{equation*}
%\left(Y^N\right)'(u_k)=\sum_{j=0}^N\xi_j\phi'_j(u_k)=\xi_k.
%\end{equation*}
%Then, the model defined in (\ref{propModel}) can be reformulated as
%\begin{equation*}
%Y^N(x)=Y^N(0)+\sum_{j=0}^N(Y^N)'(u_j)\phi_j(x),
%\end{equation*}
%since $\phi_j(0)=0, \ j=0,\ldots,N$. To ensure the pathwise uniform convergence of $Y^N$ to the original GP $Y$, see \cite{maatouk:hal-01096751}, we assume that $Y^N(0)=Y(0)$ and $(Y^N)'(u_j)=Y'(u_j), \ j=0,\ldots,N$. Then
%\begin{equation*}
%Y^N(x)=Y(0)+\sum_{j=0}^NY'(u_j)\phi_j(x)
%\end{equation*}
%is non-increasing (resp. non-decreasing) \textit{if and only if} all the coefficients $\xi_j=Y'(u_j), \ j=0,\ldots,N$ are nonpositive (resp. nonnegative). Additionally, the $(N+2)\times (N+2)$ covariance matrix $\Gamma^N$ of the Gaussian vector $\boldsymbol{\xi}$ can be written as~:
%\begin{equation*}
%\Gamma^N=\left[\begin{matrix}
%K(0,0) && \frac{\partial K}{\partial x'}(0,u_j) \\\\
%\frac{\partial K}{\partial x}(u_i,0) && \frac{\partial^2K}{\partial x\partial x'}(u_i,u_j)\\
%\end{matrix}
%\right]_{0\leq i,j\leq N},
%\end{equation*}
%where $K$ is the covariance matrix of the original GP $Y$ and $u_j, \ j=0,\ldots,N$ are the knots of the input subdivision. 
%Indeed, for all $i,j=0,\ldots,N$, we have
%\begin{eqnarray*}
%&&\cov(\xi_i,\xi_j)=\cov(Y'(u_i),Y'(u_j))=\frac{\partial^2}{\partial x\partial x'}\cov(Y(u_i),Y(u_j))=\frac{\partial^2K}{\partial x\partial x'}(u_i,u_j), \\
%&& \cov(\eta,\xi_j)=\cov(Y(0),Y'(u_j))=\frac{\partial}{\partial x'}\cov(Y(0),Y(u_j))=\frac{\partial K}{\partial x'}(0,u_j).
%\end{eqnarray*}

%\paragraph{Simulated paths.}
In order to construct curves which are compatible with market quotes, linear equality constraints like the one given in Subsection~\ref{ILC} have to be imposed on the process $Y^N$ at some points $x^{(1)},\ldots,x^{(m)}$ in $\Dom$. Then, if $Y^N(X)=\left(Y^N\left(x^{(1)}\right), \ldots,Y^N\left(x^{(m)}\right)\right)^\top$ denotes the vector of values involved in the curve construction and given Equation~\eqref{propModel}, the condition $A\cdot Y^N(X) = \boldsymbol{b}$ translates into the following linear equality constraint on the Gaussian vector $\boldsymbol{\xi}$~:
% as defined in \eqref{initLinConst} can be expressed as
%Applying the linear condition defined in Equation~(\ref{initLinConst}) to the proposed model $Y^N(x)$ defined in (\ref{propModel}), we get 
%\begin{equation}\label{EqCond}
%A\cdot Y^N(X)=A\cdot \left(\eta+\sum_{j=0}^N\xi_j\phi_j(X)\right)=\left(A\cdot \Phi\right)\cdot \boldsymbol{\xi}=\boldsymbol{b},
%\end{equation}
\begin{equation}\label{EqCond}
A \cdot \Phi \cdot \boldsymbol{\xi}=\boldsymbol{b},
\end{equation}
where $\Phi$ is a $m\times (N+2)$ matrix defined as
\begin{eqnarray*}
\Phi_{i,j}:=\left\{\begin{array}{ll}
1 & \mbox{for} \ i=1,\ldots,m \ \text{and} \ j=1,\\
\phi_{j-2}\left(x^{(i)}\right) & \mbox{for} \ i=1,\ldots,m \ \text{and} \ j=2,\ldots,N+2.
\end{array}\right.
\end{eqnarray*}
Note that, generally speaking,  the linear equality condition \eqref{EqCond} on $\boldsymbol{\xi}$ admits solutions only when $N+2\geq n$ as  $A \cdot \Phi$ is a matrix of dimension $n\times (N+2)$.

\paragraph{Curve simulation.}  Conditional GP satisfying both monotonicity and linear equality constraints can be sampled by  generating truncated Gaussian vector $\boldsymbol{\xi}$ restricted to~: 
\begin{eqnarray*}\left\{\begin{array}{ll}
B \cdot \boldsymbol{\xi}=\boldsymbol{b} & \qquad \mbox{linear equality condition}\\
\boldsymbol{\xi}\in \Cxi & \qquad \mbox{monotonicity constraint}
\end{array}\right.
\end{eqnarray*}  
where,  $B  = A \cdot \Phi$ and for instance $\Cxi=\left\{\boldsymbol{\xi}\in\mathbb{R}^{N+2} \ : \ \xi_j\leq 0, \ j=0,\ldots,N\right\}$ for non-increasing constraints. 
 Then, simulated paths can be sampled in two steps.  First, the conditional  distribution of the vector $\boldsymbol{\xi}$ given $B \cdot \boldsymbol{\xi}=\boldsymbol{b}$ is still Gaussian with mean 
 $$
 (B\Gamma^N)^\top\left(B\Gamma^N  B^\top\right)^{-1}\boldsymbol{b}
 $$
 and covariance matrix
 $$
 \Gamma^N-\left(B\Gamma^N\right)^\top\left(B\Gamma^N  B^\top\right)^{-1}B\Gamma^N,
 $$
so that it can be simulated very efficiently.
% \hassansuppr{\did{, using the same kind of results as in Section~\ref{ILC}}}.
%\hassan{\begin{equation}\label{IntXi}
%\xi \suchthat (A\cdot\Phi)\cdot\xi=\boldsymbol{b} \sim \mathcal{N}\left((A\Phi\Gamma^N)^\top\left(A\Phi\Gamma^N\Phi^\top A^\top\right)^{-1}\boldsymbol{b},\Gamma^N-(A\Phi\Gamma^N)^\top\left(A\Phi\Gamma^N\Phi^\top A^\top\right)^{-1}A\Phi\Gamma^N\right).
%\end{equation}
Then, the simulation of truncated Gaussian vectors restricted to, for instance, negativity of the components (here $\boldsymbol{\xi}\in \Cxi$) can be done by using improved rejection sampling algorithm such as the one described in \cite{MaatoukMCQMC2014} and \cite{Robert}.
%Then Gaussian vectir  $\boldsymbol{\xi}$ with the above conditional distribution~\eqref{IntXi}}. Finally, by using improved rejection sampling \cite{MaatoukMCQMC2014}, we select the non-positive coefficients ($\xi_j\leq 0$) for non-increasing constraints}.
 By Equation~\eqref{propModel} we get sample paths that fulfil both constraints. These simulations can be used to construct confidence intervals for the value of the curve at each point $x$ in $\Dom$.

\paragraph{Most likely curve.} 

Given a covariance kernel $K$ and its estimated parameters, it is possible to determine the most likely path of the conditional Gaussian process under both linear and monotonicity constraints. This most likely curve corresponds to the mode\footnote{The maximum of the probability density function.} of the finite-dimensional truncated Gaussian vector $\boldsymbol{\xi}$ \citep[see also][for a discussion about the mode]{abrahamsen2001kriging}.
%This estimator corresponds to a particular spline (see~\cite{bay:hal-01136466}). 
%Using the model defined in~(\ref{propModel}), \did{this} new estimator called mode is derived. 
In bayesian statistics, it is known as the Maximum A Posteriori (MAP) estimator \citep[see][for more details]{maatouk:hal-01096751}. Its  expression is given by~:
\begin{equation}\label{propEstimator}
M^N_{K}\left(x \suchthat A,\boldsymbol{b}\right)=\nu+\sum_{j=0}^N\nu_j\phi_j(x),
\end{equation}
where $\boldsymbol{\nu}=(\nu,\nu_0,\ldots,\nu_N)^\top\in\mathbb{R}^{N+2}$ is the solution of the following convex optimization problem~:
\begin{equation}\label{eqNU}
\boldsymbol{\nu}=\arg\min_{\boldsymbol{c}\in \Cxi\cap \Ixi(A,\boldsymbol{b})}\left(\frac{1}{2}\boldsymbol{c}^\top\left(\Gamma^N\right)^{-1}\boldsymbol{c}\right),
\end{equation}
and where $\Gamma^N$ is the covariance matrix of the Gaussian vector $\boldsymbol{\xi}$ defined in \eqref{eq:cov_mat_one_dim}. The vector $\boldsymbol{\nu}$ can be seen as the mode of the Gaussian vector $\boldsymbol{\xi}$ restricted to $\Cxi\cap \Ixi(A,\boldsymbol{b})$, where $\Cxi$ is the set of vectors satisfying  monotonicity constraints and $\Ixi(A,\boldsymbol{b}) = \left\{\boldsymbol{\xi}\in\mathbb{R}^{N+2} \ : \ A\cdot \Phi\cdot\boldsymbol{\xi}=\boldsymbol{b}\right\}$ is the set of vectors which are compatible with the linear equality constraint.
%$ represents the \hassan{linear equality condition} $A\cdot \Phi\cdot\boldsymbol{\xi}=\boldsymbol{b}$. 
Obviously, this curve satisfies  both constraints. 
%\hassan{linear equality condition}~\eqref{EqCond} and the monotonicity constraints in the entire domain. 
Additionally, it does not depend on the variance  $\sigma^2$ of the  covariance function $K$ since $\sigma^2$ is a multiplicative constant in the matrix $\Gamma^N$ and then does not affect the $\arg\min$ in Equation~\eqref{eqNU}. In \cite{bay:hal-01136466} and \cite{2016arXiv160202714B}, the convergence of  the proposed estimator~(\ref{propEstimator})  as $N$ tends to infinity is studied and its limit corresponds to a constrained spline function (which depends on the underlying kernel function). By this methodology, one can thus retrieve classical spline interpolation with the additional possibility of getting confidence intervals.

\subsection{Two dimensional case}\label{TDC}

%The previous construction method can be easily extended 
As explained in Remark~\ref{rem:surface}, the term-structure construction can be stated in a two-dimensional setting. The aim is to incorporate in the curve construction process  market information at several quotation dates. Contrary to the previous section, the output is now a surface which may represent the evolution of a reference quantity with respect to time-to-maturities and quotation dates. More formally, we consider a two-dimensional input variable $\boldsymbol{x}=(x,t)$ which is assumed to belong to a rectangle $\Dom = [\underline{x}, \overline{x}]\times [\underline{t}, \overline{t}]$ of $\R^2$. The variable $x$ may represent time-to-maturities whereas the variable $t$ may represent evolution of time or quotation dates. For the sake of simplicity, we assume that the unknown bivariate real function $f$ is monotone, say non-increasing, with respect to the first input variable only~:
%\begin{equation}
%\label{eq:bivariate_monotone_constraint}
%x_{1}\leq x_{2} \quad \Rightarrow \quad f(x_{2},t)\leq f(x_{1},t), \qquad  \mbox{ for all }  t \in[\underline{t}, \overline{t}],\;  x_{1}, x_{2} \in[\underline{x}, \overline{x}].  
%\end{equation}
\begin{equation}
\label{eq:bivariate_monotone_constraint}
x_{a}\leq x_{b} \quad \Rightarrow \quad f(x_{b},t)\leq f(x_{a},t), \qquad  \mbox{ for all }  t \in[\underline{t}, \overline{t}],\;  x_{a}, x_{b} \in[\underline{x}, \overline{x}].  
\end{equation}

More general inequality conditions can be considered in this bivariate setting \citep[see][for more details]{maatouk:hal-01096751}. The aim of this section is to explain how to construct a process that simultaneously  satisfies a series of linear equality constraints (one for each considered quotation date) and a monotonicity constraint as described in~\eqref{eq:bivariate_monotone_constraint}. As in the one-dimensional setting, we start with an original bivariate Gaussian process $Y$ with zero-mean and with a covariance function $K$. 
%and we assume for simplicity, we assume that $Y$ is a zero-mean  GP. 
\\

The idea is the same as the one dimensional case presented in Section~\ref{ODC}. We begin by discretizing the input rectangle $\Dom$ in a $(N_x+1)\times (N_t+1)$ grid which for simplicity is assumed to be regular. 
%The interval  a regular subdivision $u_0<\ldots < u_N$ with $u_0=\underline{x}$, $u_N=\overline{x}$ and with a constant mesh $\delta$, so that $u_j=u_0+j\delta$, $j=0, \ldots, N$.
The subdivision of the $x$-axis is $u_0<\ldots < u_{N_x}$ with $u_0=\underline{x}$, $u_{N_{x}}=\overline{x}$ and a constant mesh $\delta_{x}$, so that $u_i=u_0+i\delta_{x}$, $i=0, \ldots, N_{x}$. The subdivision of the $y$-axis is $v_0<\ldots < v_{N_t}$ with $v_0=\underline{t}$, $v_{N_{t}}=\overline{t}$ and a constant mesh $\delta_{t}$, so that $v_j=v_0+j\delta_{t}$, $j=0, \ldots, N_{t}$.
%where $u_i=u_0+i\delta_x$, $i=0, \ldots, N_x$ and the   subdivision of the $t$-axis is $v_0<\ldots < v_{N_t}$ where $v_j=v_0+j\delta_t$, $j=0, \ldots, N_t$, with respective \ares{constant} mesh sizes $\delta_x$, $\delta_t$ $\in \R_+^*$. 
The following developments  can be easily extended to irregular grids. As in the one-dimensional case, we consider an original Gaussian  process $Y$ with covariance function $K$. We then  define the finite-dimensional approximation of $Y$ as the process $Y^N$ such that 
\begin{equation}\label{2Dmodel}
Y^N(x,t)=\sum_{i=0}^{\Nx}\sum_{j=0}^{\Nt}\xi_{i,j}g_i(x)h_j(t), \qquad  \mbox{ for all }  (x,t)\in\Dom,
\end{equation} 
where $g_i(x)=\max\left(1-\frac{|x-u_i|}{\delta_x},0\right)$ and $h_j(t)=\max\left(1-\frac{|t-v_j|}{\delta_t},0\right)$ are  hat functions centered at the knots $u_i$ and $v_j$, for $i=0, \ldots, N_x$ and $j=0, \ldots, N_t$ and $\boldsymbol{\xi}=(\xi_{0,0},\xi_{0,1},\ldots,\xi_{i,j},\ldots,\xi_{\Nx,\Nt})^\top$ is a zero-mean Gaussian vector with components $(\boldsymbol{\xi})_{\rho_{ij}}=\xi_{i,j}$ where $\rho_{ij}=(\Nt+1)i+j+1$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$.  Let $N_{tot}=(N_x+1) (N_t+1)$ be the size of the column vector $\boldsymbol{\xi}$. Choosing  $\xi_{i,j} = Y(u_i,v_j)$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$ guarantees that the finite-dimensional  process $Y^N$ converges on $\Dom$ towards $Y$ almost surely as $N_x$ and $N_t$ tends to infinity \citep[see][]{maatouk:hal-01096751}. In that case, the covariance matrix $\Gamma^N\in\mathbb{R}^{N_{tot}^2}$ of the Gaussian vector $\boldsymbol{\xi}$ can be written as~:
\begin{equation*}
\Gamma^N_{\rho_{ij},\rho_{i'j'}}=\cov(\xi_{i,j},\xi_{i',j'})= K\left((u_i,v_j),(u_{i'},v_{j'})\right),
\end{equation*}
where $i,i'=0,\ldots,\Nx$ and $j,j'=0,\ldots,\Nt$. As can be shown in the next proposition, the monotonicity constraint~\eqref{eq:bivariate_monotone_constraint} reduces to a linear inequality condition on the vector $\boldsymbol{\xi}$. 
\begin{pro}[Monotonicity]
\label{pro:monotone_2D}
The process $Y^N$ defined in Equation~\eqref{2Dmodel} is non-increasing (resp. non-decreasing) with respect to the first variable $x$ \textit{if and only if} all the coefficients $\xi_j, \ j=0,\ldots,N$ are such that
\begin{equation}\label{Dim2LinConst}
\xi_{i-1,j}\leq  \xi_{i,j},\; (\text{resp. } \xi_{i-1,j}\geq  \xi_{i,j}) \qquad i=1,\ldots,\Nx \ \mbox{and} \ j=0,\ldots,\Nt.
\end{equation}
\end{pro}
\begin{proof}
The proof is similar to the one of Proposition~\ref{pro:monotone_2D}.
\end{proof}

%The monotonicity (non-increasing) constraints on $Y^N$ with respect to the first input variable are equivalent to the following linear inequality constraints on the coefficients~:
%\begin{equation}\label{Dim2LinConst}
%\xi_{i-1,j}\leq \xi_{i,j}, \qquad i=1,\ldots,\Nx \ \mbox{and} \ j=0,\ldots,\Nt.
%\end{equation}


%When $\xi_{i,j} = Y(u_i,v_j)$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$, the uniform convergence of $Y^N$ towards $Y$ is stated in \cite{maatouk:hal-01096751}. It can also be shown that $Y^N$ is again a Gaussian process and that its covariance function $K_N$ uniformly converges to the covariance function $K$ of $Y$.

%\paragraph{Convergence.} In that case, the basis functions verify~:
%\begin{equation*}
%g_i(u_k)h_j(v_l)=\indic{i=k}\indic{j=l}.
%\end{equation*}
%As a consequence, $Y^N(u_k,v_l)=\xi_{k,l}$. One thus can choose $\xi_{i,j}$ such that $Y^N(u_k,v_l)=Y(u_k,v_l), \ k=0,\ldots,\Nx, \ l=0,\ldots,\Nt$ and the model defined in (\ref{2Dmodel}) can be reformulated as
%\begin{equation*}
%Y^N(x,t)=\sum_{i=0}^{\Nx}\sum_{j=0}^{\Nt}\xi_{i,j}g_i(x)h_j(t)=\sum_{i=0}^{\Nx}\sum_{j=0}^{\Nt}Y(u_i,v_j)g_i(x)h_j(t),
%\end{equation*} 
%which converges uniformly to the original GP $Y$, see \cite{maatouk:hal-01096751}.
%In order to express the linear constraints on $Y^N$, we define the column vector $\boldsymbol{\xi}=(\xi_{0,0},\xi_{0,1},\ldots,\xi_{i,j},\ldots,\xi_{\Nx,\Nt})^\top$ with components $(\boldsymbol{\xi})_{\rho_{ij}}=\xi_{i,j}$ where $\rho_{ij}=(\Nt+1)i+j+1$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$. The covariance matrix $\Gamma^N\in\mathbb{R}^{\left((\Nx+1)\cdot (\Nt+1)\right)^2}$ of this vector $\boldsymbol{\xi}$ can be written as~:
%\begin{equation*}
%\Gamma^N_{\rho_{ij},\rho_{i'j'}}=\cov(\xi_{i,j},\xi_{i',j'})=\cov\left(Y(u_i,v_j),Y(u_{i'},v_{j'})\right)=K\left((u_i,v_j),(u_{i'},v_{j'})\right),
%\end{equation*}
%where $i,i'=0,\ldots,\Nx$ and $j,j'=0,\ldots,\Nt$.


In order to construct surfaces which are compatible with market quotes observed at times  $t=t_1, \ldots, t_{I}$, $I$ different linear equality constraints like the one given in Section~\ref{ILC} has to be imposed on the bivariate process $Y^N$. Let us consider that the market-fit conditions involve $m$ points $x^{(1)},\ldots,x^{(m)}$ for any quotation time $t=t_1, \ldots, t_I$. Then, if $Y^N(X,t)$ $=$ $\left(Y^N\left(x^{(1)},t\right), \ldots,\right.$ $\left.Y^N\left(x^{(m)},t\right)\right)^\top$ denotes the vector of values involves in the curve construction at time $t$ and given Equation~\eqref{2Dmodel}, the condition $A_t\cdot Y^N(X,t)=\boldsymbol{b}_t$ translates into the following linear equality constraint on the Gaussian vector $\boldsymbol{\xi}$~:
\begin{equation}
\label{EqCond_2D}
A_t\cdot H_t\cdot \boldsymbol{\xi}=\boldsymbol{b}_t
\end{equation}
where the $m\times N_{tot}$ matrix $H_t$ has components $(H_t)_{k,\rho_{ij}}=g_i\left(x^{(k)}\right)h_j(t)$, $k=1, \ldots, m$ and $\rho_{ij}=(\Nt+1)i+j+1$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$. The former condition has to hold simultaneously for every time $t=t_1, \ldots, t_I$, which can be summarized as one single linear equality constraint  
\begin{equation}
\label{EqCond_2D_bis}
B \cdot \boldsymbol{\xi} = \boldsymbol{b},
\end{equation}
where the $nI \times N_{tot}$ matrix $B$ (resp.  $\boldsymbol{b}$) is formed by  vertical concatenation of matrices $A_t \cdot H_t$ (resp. $\boldsymbol{b}_t$) for $t=t_1, \ldots, t_I$.
%The matrix $H_t$ has size $m\times N_{tot}$ and the column vector $\boldsymbol{\xi}$ has size $N_{tot}$, where $N_{tot}=(N_x+1) (N_t+1)$. 
Notice that, generally speaking, this linear system  admits  solutions  when $N_{tot} \geq nI$.
%when $N_{x} +1 \geq m$ and $N_{t} +1 \geq T$.

%We moreover assume that function $f$ satisfies a set of  linear equality constraints given as follows~: 
%\begin{equation*}
%A_t\cdot Y^N(X,t)=\boldsymbol{b}_t, \qquad \ t=t^{(1)}, \ldots,t^{(T)},
%\end{equation*}
%where $Y^N(X,t)=\left(Y^N(x^{(1)},t),\ldots, Y^N(x^{(m)},t) \right)$.  


%The monotonicity (non-increasing) constraints on $Y^N$ with respect to the first input variable are equivalent to the following linear inequality constraints on the coefficients~:
%\begin{equation}\label{Dim2LinConst}
%\xi_{i-1,j}\leq \xi_{i,j}, \qquad i=1,\ldots,\Nx \ \mbox{and} \ j=0,\ldots,\Nt.
%\end{equation}


\paragraph{Surface simulation.} 
As in the one-dimensional setting, the simulation of the bivariate Gaussian process $Y^N$ conditionally to both linear equality  and monotonicity constraints reduces to the simulation of the Gaussian vector $\boldsymbol{\xi}$ restricted to 
\begin{eqnarray*}\left\{\begin{array}{ll}
B\cdot \boldsymbol{\xi}=\boldsymbol{b} & \qquad \mbox{linear equality condition}\\
\boldsymbol{\xi}\in \Cxi & \qquad \mbox{monotonicity constraint}
\end{array}\right.
\end{eqnarray*}
%\begin{eqnarray*}\left\{\begin{array}{ll}
%B_t\cdot \boldsymbol{\xi}=\boldsymbol{b}_t, \quad \mbox{for all }t=t^{(1)}, \ldots,t^{(T)} & \qquad \mbox{\hassan{linear equality conditions}}\\
%\boldsymbol{\xi}\in C_{\xi} & \qquad \mbox{monotonicity constraints}
%\end{array}\right.
%\end{eqnarray*}
where, for instance, $\Cxi=\left\{\boldsymbol{\xi}\in\mathbb{R}^{N_{tot}} \ : \ \xi_{i-1,j}\leq \xi_{i,j}\right\}$ for a non-increasing constraint. 
%The matrix $H_t$ has components $(H_t)_{k,\rho_{ij}}=g_i(x_k)h_j(t)$, where $\rho_{ij}=(\Nt+1)i+j+1$, $i=0,\ldots,N_x$ and $j=0,\ldots,N_t$. 
%The matrix $H_t$ has size $n\times N_{tot}$ and the column vector $\boldsymbol{\xi}$ has size $N_{tot}$, where $N_{tot}=(N_x+1) (N_t+1)$. Notice that, generally speaking, all the linearity conditions admit here solutions only when $N_{x} +1 \geq m$ and $N_{t} +1 \geq T$.\\
The simulation procedure is the same as the one presented in Section~\ref{ODC}. 
These simulations can be used to construct confidence intervals for the surface value at each point $(x,t)$ in $\Dom$.

%, using the fact that all constraints $\left(A_t\cdot H_t\right)\cdot \boldsymbol{\xi}=\boldsymbol{b}_t$, $t=t^{(1)}, \ldots,t^{(T)}$ can be summarized as one single matrix constraints on the column vector $\boldsymbol{\xi}$.

%$(T\cdot n) \times ((\Nx+1)\cdot (\Nt+1))$

\paragraph{Most likely surface.} The methodology can be easily extended from Section~\ref{ODC}.

\subsection{Parameters estimation}\label{PE}

The most likely path of the constrained process (mode estimator) 
depends  on the choice of the underlying Gaussian process or equivalently on its covariance function $K$. In this section, we investigate the estimation of its parameter, i.e., the length and the variance hyper parameters. 
%We now investigate the estimation of covariance hyper parameters of Gaussian processes with monotonicity constraints is considered. 
In the literature,  estimation of the covariance function hyper parameters is usually done using two types of methods~: Maximum Likelihood (ML) estimators as in \cite{sant:will:notz:2003} and Cross Validation (CV) methods  as in \cite{Bachoc201355}, \cite{Cressie1993Wiley} and \cite{Roustant:Ginsbourger:Deville:2012:JSSOBK:v51i01}. Both methods ML and CV are not suited to monotonicity constraints.\\

Recently, an Adapted Cross-Validation (ACV) technique has be proposed by  \cite{Maatouk201538} to estimate covariance hyper-parameters of Gaussian processes in the presence of inequality constraints. The main idea is to consider the mode curve (as opposed to the mean curve) as the estimator to be used in the cross-validation method. 
%In the following, we focus on these ACV methods for estimating length and variance hyper-parameters.
The principle of cross-validation is to select the set of parameters that minimizes a distance between observed values and their estimates while successively omitting some set of observations. 
Thus, cross-validation is similar in spirit to backtesting \citep[see][]{kerkhof2004backtesting}, but omitted data are not necessarily taken in chronological order.


\paragraph{Length parameters.}

In classical kriging, the usual cross-validation estimator of the covariance  length parameter $\boldsymbol{\theta}$  is constructed from the so-called Leave One Out (LOO) mean square error criterion. Given that the unknown function $f$ takes values $y_1, \ldots, y_n$ at points $x^{(1)}, \ldots, x^{(n)}$ (pure interpolation constraints), the cross-validation estimator $\hat{\boldsymbol{\theta}}_{CV}$ of  $\boldsymbol{\theta}$ is defined as 
\begin{equation}
\hat{\boldsymbol{\theta}}_{CV}=\arg\min_{\boldsymbol{\theta}\in\Theta}\sum_{i=1}^n\left(y_i-\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})\right)^2,
\end{equation}
where $\boldsymbol{y}_{-i}=(y_1,\ldots,y_{i-1},y_{i+1},\ldots,y_n)^\top$ and $\Theta$ is a compact subset of $\mathbb{R}^d$. The estimator 
$\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})=\mathbb{E}\left(Y(x^{(i)})\suchthat Y(X^{(-i)})=\boldsymbol{y}_{-i}\right)$ is the kriging mean at point  $x^{(i)}$ obtained by removing observation $y_i$ in the estimation process (see Equation \eqref{eq:kriging_mean}). The vector $Y\left(X^{(-i)}\right)$ is the same vector as $Y(X)$ without component $Y\left(x^{(i)}\right)$.\\




%The estimate of the length hyper parameters proposed in \cite{kerkhof2004backtesting} is an extension of the classical 
%
%In this paragraph, our aim is to estimate the length hyper parameters $\boldsymbol{\theta}$ of the covariance function $K$. Let us recall that the Leave One Out (LOO) mean square error criterion is defined as an estimator of the correlation length hyper parameter $\boldsymbol{\theta}$~:
%\begin{equation*}
%\hat{\boldsymbol{\theta}}_{CV}=\arg\min_{\boldsymbol{\theta}\in\Theta}\sum_{i=1}^n\left(y_i-\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})\right)^2,
%\end{equation*}
%where $\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})$ is the kriging mean \eqref{eq:kriging_mean} without the point $(x^{(i)}, y_i)$.


%with \did{$\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})=\mathbb{E}\left(Y(x)\suchthat Y(x^{(1)})=y_1,\ldots, Y(x^{(i-1)})=y_{i-1},Y(x^{(i+1)})=y_{i+1},\ldots,Y(x^{(n)})=y_n\right)$,}


%with \did{$\hat{y}_{i,\boldsymbol{\theta}}(\boldsymbol{y}_{-i})=\mathbb{E}\left(Y(x)\middle | Y\left(x^{(1)}\right)=y_1,\ldots, Y\left(x^{(i-1)}\right)=y_{i-1},Y\left(x^{(i+1)}\right)=y_{i+1},\ldots,Y\left(x^{(n)}\right)=y_n\right)$,
% \did{$\boldsymbol{y}_{-i}=(y_1,\ldots,y_{i-1},y_{i+1},\ldots,y_n)^\top$} and $\Theta$ is a compact subset of $\mathbb{R}^d$. At any location $x_i$, the criterion compares the given real values $y_i$ with their estimation given all other observations. 
%
 
The previous criterion cannot be used in the presence of monotonicity constraints since 
the classical kriging mean estimate does not respect such kind of constraints.  As suggested in \cite{Maatouk201538} under pure interpolation constraints,  the most likely curve (mode estimator) defined in Section \ref{ODC} can be used instead of the kriging mean since the former satisfies the monotonicity constraints. Then, a new LOO criterion adapted to monotonicity constraints can be defined as~:
\begin{equation}\label{LOO}
\hat{\boldsymbol{\theta}}_{ACV}=\arg\min_{\boldsymbol{\theta}\in\Theta}\sum_{i=1}^n\left(y_i-M^N_{K}\left(x^{(i)} \suchthat I_{n-1},\boldsymbol{y}_{-i}\right)\right)^2,
\end{equation} 
where $M^N_{K}\left(x^{(i)} \suchthat I_{n-1},\boldsymbol{y}_{-i}\right)$ is the mode estimator (defined in Equation~(\ref{propEstimator})) of $y_i$ based on all observations but $y_i$. 
As explained in  Section \ref{ODC}, the latter does not depend on the variance $\sigma^2$. The matrix
%vector $\boldsymbol{y}_{-i}=(y_1,\ldots,y_{i-1},y_{i+1},\ldots,y_n)^\top$ and 
$I_{n-1}$ is the identity matrix with dimension $n-1$ (the same as the vertical dimension as the column vector $\boldsymbol{y}_{-i}$). This methodology cannot be applied immediately in our setting since the output values $y_i$ in the LOO criterion (\ref{LOO})
are not available under linear equality constraints as defined in (\ref{EqCond}). To this end, we consider the following formulation of the LOO criterion~:
\begin{equation}\label{LOOLConst}
\hat{\boldsymbol{\theta}}_{ACV}=\arg\min_{\boldsymbol{\theta}\in\Theta}\sum_{i=1}^{n}\left(b_i-\left(A\cdot M^N_{K}\left(X \suchthat A_{-i},\boldsymbol{b}_{-i}\right)\right)_i\right)^2,
\end{equation}
where the subscript $i$ refers to the $i$-th component, $A_{-i}$ and $\boldsymbol{b}_{-i}$ are respectively the matrix and the vector without the $i^{\text{th}}$ row.\\


\paragraph{Variance parameter.}
%Applying the same kind of cross validation techniques to the variance hyper parameter $\sigma$ (see~\cite{Bachoc201355} and \cite{maatouk:hal-01153960} for details), we get

In the case of pure interpolation constraints and without monotonicity constraints, several classical criterions can be used for the estimation of $\sigma$. For instance, in \cite{Bachoc201355}, an estimator $\hat{\sigma}_{CV}$ is defined as the parameter $\sigma$ such that the following equality holds 
%\begin{equation}\label{sigSCV}
%\hat{\sigma}_{CV}^2=\frac{1}{n}\sum_{i=1}^n\frac{\left(y_i-\hat{y}_{i,\hat{\boldsymbol{\theta}}_{CV}}(\boldsymbol{y}_{-i})\right)^2}{\hat{c}^2_{i,-i}},
%\end{equation}
%where $\hat{c}^2_{i,-i}$ is the kriging variance without the $i^{\text{th}}$ observation data, i.e., $\var\left(Y(x^{(i)})\suchthat Y(\boldsymbol{x}^{(-i)})=\boldsymbol{y}_{-i}\right)$, where $\sigma$ if fixed to 1.
\begin{equation}\label{sigSCV}
\frac{1}{n}\sum_{i=1}^n\frac{\left(y_i-\hat{y}_{i,\hat{\boldsymbol{\theta}}_{CV}}(\boldsymbol{y}_{-i})\right)^2}{\mathds{E}\left(\left(Y\left(x^{(i)}\right)-\hat{y}_{i,\hat{\boldsymbol{\theta}}_{CV}}(\boldsymbol{y}_{-i})\right)^2\suchthat Y(X^{(-i)})=\boldsymbol{y}_{-i}\right)}=1.
\end{equation}
Given our linear equality conditions and monotonicity constraints, we instead propose to define $\hat{\sigma}_{ACV}$ as the parameter $\sigma$ such that the following equality holds
%\begin{equation}\label{sigSCV_2}
%\frac{1}{n}\sum_{i=1}^n\frac{\left(b_i-\left(A\cdot M^N_{K}\left(X \suchthat A_{-i},\boldsymbol{b}_{-i}\right)\right)_i\right)^2}{\frac{1}{n_{sim}}\sum_{j=1}^{n_{sim}}\left(AY_{j,-i}\left(X\right)-A M^N_{K}\left(X\suchthat A_{-i},\boldsymbol{b}_{-i}\right)\right)_i^2}=1,
%\end{equation}
%where $Y_{j,-i}\left(X^{(i)}\right)$ represents the value of a $j^{\text{th}}$ simulation without the $i^{\text{th}}$ observation at $x^{(i)}$ and $n_{sim}$ represents the number of simulation.
\begin{equation}\label{sigSCV_2}
\frac{1}{n}\sum_{i=1}^n\frac{\left(b_i-\left(A\cdot M^N_{K}\left(X \suchthat A_{-i},\boldsymbol{b}_{-i}\right)\right)_i\right)^2}{\mathds{E}\left(\left(AY\left(X\right)-A M^N_{K}\left(X\suchthat A_{-i},\boldsymbol{b}_{-i}\right)\right)_i^2 \suchthat \mathcal{D}_i \right)}=1,
\end{equation}
where $\mathcal{D}_i$ is the set of monotonicity constraints on the whole domain and the additional linear constraints without the $i^{\text{th}}$ one, i.e. $A_{-i} Y(X^{(-i)}) = \boldsymbol{b}_{-i}$ where $A_{-i}$ is the matrix $A$ without the $i$-th row.
The expectation is estimated by simulation, approximating the process $Y$ by $Y^N$.

%where $Y_{j,-i}\left(X^{(i)}\right)$ represents the value of a $j^{\text{th}}$ simulation without the $i^{\text{th}}$ observation at $x^{(i)}$ and $n_{sim}$ represents the number of simulation.

%\textsc{\begin{equation}\label{eq:estim_y}
%\hat{y}_{-i}(\hat{\boldsymbol{\theta}}_{ACV}, \sigma) = \mathbb{E} \left(Y(X)\suchthat A_{-i}Y(X^{(-i)})=\boldsymbol{b}_{-i}\right)
%\end{equation}
%and where 
%\begin{equation}\label{eq:estim_c}
%\hat{c}^2_{i,-i} (\hat{\boldsymbol{\theta}}_{ACV}, \sigma)= \mathbb{E} \left(\left(\left[A Y(X)\right]_i -\left[A\hat{y}_{-i}(\hat{\boldsymbol{\theta}}_{ACV}, \sigma) \right]_i\right)^2 \suchthat A_{-i}Y(X^{(-i)})=\boldsymbol{b}_{-i}\right)
%\end{equation}
%is the variance of $\left[A Y(X)\right]_i$ given $A_{-i}Y(X^{(-i)})=\boldsymbol{b}_{-i}$, where the hyper-parameters of the covariance function of $Y$ are $\hat{\boldsymbol{\theta}}_{ACV}$ and  $\sigma$.
%}The idea is the same for estimating the variance hyper parameter in Equation~(\ref{sigSCV}).

\section{Empirical investigation}\label{EmpInv}

In this section, the construction method developed in Section \ref{IntModIConst} is illustrated in different financial applications. Based on  market quotes observed at different  dates, we construct curves together with confidence intervals for term-structures of OIS discount rates, term-structure of zero-coupon swaps rates and term-structure of CDS implied default probabilities. Using the bivariate setting of Section \ref{TDC}, we also build interest-rate and default probability surfaces by considering time (quotation dates) as an additional dimension.


%on real data
% given by the linear condition defined in (\ref{EqCond}) to estimate term-structure functions. 
%As explain in Section \ref{Sec:term-structure} provides an incomplete information on the curve in the form of a market-fit condition.
% In general, the interpolation conditions or given data are seen as design points. This is not the case in the estimating of discount factors, where the real function at some particular points is known to satisfy a system of equations.

%to pass through some intervals see ARESKI2014. It is interpolated at some design points in the first part of the input domain and then it respects some linear conditions in the second part. In this framework, the input set represents the years and the output corresponds to the values of discount factors.

\subsection{Interest-rate curves based on Swaps versus Euribor}\label{CD}%\paragraph{The data.} 


We apply the kriging procedure to construct zero-coupon swap curves  based on market quotes of fixed-vs-floating interest-rate swaps   for different standard maturities.
Market observations are given as par rates of Swaps vs Euribor 6M. We consider the 10 quotation dates given in Table~\ref{EPOV}. For each quotation date, the term-structure is built  on $14$ swap rates $S_1, \ldots, S_{14}$ associated with standard maturities in years belonging to the set $E:=\{1,\ldots,10,15,20,30,40\}^\top$. As explained in Section~\ref{Sec:term-structure}, each observed par rate provides an indirect  information on the curve. This information takes the
  form of a linear relation given by~\eqref{eq:OIS}. For each standard maturity $T\in E$, this relation involves the value $P(t,k)$ of the curve at time horizon $k=1, \ldots, T$.
As a result, the curves are compatible with observed quotes if, for each observation date $t$, the vector of discount factors $P(t, X) := (P(t, 1), \ldots, P(t,40))^{\top}$ satisfies a linear system of the form 
\begin{equation}
\label{eq:market_fit_Swaps_Euribor}
A_t \cdot P(t, X)=\boldsymbol{b}_t,
\end{equation}
where $A_t$ is  a $14\times 40$ real matrix and $\boldsymbol{b}_t= (1,\ldots,1)^\top\in \mathbb{R}^{14}$. In this case, we have $n=14$ products whose value depends on $m=40$ points of the curve.\\

We consider that the associated discount factor curve $P(t,X)$ is one realization of a decreasing spatial process which starts from $1$ ($P(t,0)=1$) and satisfies the linear condition~\eqref{eq:market_fit_Swaps_Euribor}.  Section~\ref{IntModIConst} explains how to construct and simulate a  process with monotonicity and linear equality constraints. This construction involves a finite-dimensional approximation $Y^N$ of Gaussian processes as defined in~\eqref{propModel}.  The latter depends on a $N+2$-dimensional Gaussian vector $\boldsymbol{\xi}$ and  a set of basis functions $\phi$ defined on a subdivision of the input domain $D=[0,40]$. We consider here a regular subdivision $u_j=u_0+j\delta$, $j=0, \ldots, N$, where $u_0=0$, $N=50$ sub-intervals and $\delta=1$.
%\ares{Ares : Describe the  discretisation of the input interval $\Dom$ using a finite-subdivision $u_0<\ldots < u_N$ where $u_j=u_0+j\delta$, $j=0, \ldots, N$ and $\delta$ is the mesh of this subdivision. }
Since the curve $T\rightarrow P(t,T)$ is known to start from 1, i.e., $P(t,0)=1$, the linear equality condition on $\boldsymbol{\xi}$ defined in (\ref{EqCond}) can be reformulated as follows~:
\begin{equation}\label{EqXiCond}
\left(\begin{matrix}
1 & \phi_0(0) & \ldots & \phi_N(0)\\
& & A_t\cdot \Phi& &  
\end{matrix}
\right)\boldsymbol{\xi}=\left(\begin{matrix}
1\\
\boldsymbol{b}_t
\end{matrix}\right).
\end{equation}
Then the simulation of the GP $Y^N$ conditionally to the linear equality condition and non-increasing constraints is equivalent to generating a truncated zero-mean Gaussian vector $\boldsymbol{\xi}$ restricted to Equation~(\ref{EqXiCond}) and to non-positive components.
%$C_{\xi}$, where $C_{\xi}=\left\{\boldsymbol{\xi}^\top\in\mathbb{R}^{N+2} \ : \ \xi_j\leq 0, \ j=0,\ldots,N\right\}$. 



\paragraph{Parameters estimation.}
We first consider the correlation-length hyper parameter $\theta$. It has been estimated for each quotation dates in Table~\ref{EPOV}  by using the adapted cross-validation (ACV) method described in Section \ref{PE}. We consider covariance functions of the form $K(x,x')=\sigma^2 C(x-x',\theta)$ and we discuss two alternative kernels $C$, that is the Gaussian and the Mat\'ern 5/2 kernel (see Table~\ref{kernel}).
We denote  by $\hat{\theta}_G$ (resp. $\hat{\theta}_M$) the estimated length parameter of the Gaussian (resp. Mat\'ern 5/2) covariance function. 
As shown in Table~\ref{EPOV}, the estimated length parameter $\hat{\theta}$ remains stable for both Gaussian and Mat\'ern 5/2 covariance functions (value around $25$ for $\hat{\theta}_G$ and $30$ for $\hat{\theta}_M$). Note also that the minimal value of the objective function in the LOO criterion \eqref{LOOLConst} is slightly smaller when using a Gaussian covariance function. However, as can be seen in Figure~\ref{LOOMatGauss}, the objective function (cross-validation error)  reach $0.12$ for the Gaussian covariance function whereas it never goes above $0.02$ for the Mat\'ern 5/2 covariance function. Additionally, the global minimum is easier to find in the Mat\'ern 5/2 case.

\begin{table}[hptb]
\centering
\caption{Parameters estimation using ACV methods (Swap versus Euribor 6M).}
\label{EPOV}
\begin{tabular}{ccccc}
\hline \hline
Date  &  $\hat{\theta}_G$ &  $\hat{\theta}_M$ & Gaussian optimal value & Mat\'ern 5/2 optimal value  \\   
\hline
02/06/2010   &  25.8   &  30.8  & 4.0e-06   &  \textbf{1.1e-06}   \\
05/07/2010   &  26.2   &  25.2  & \textbf{3.8e-06}   &  1.1e-05   \\
03/08/2010   &  26.6   &  27.0  & \textbf{4.1e-06}   &  4.7e-06   \\
29/11/2010   &  23.5   &  39.9  & \textbf{5.4e-07}   &  2.3e-05   \\ 
30/12/2010   &  27.5   &  28.0  & 4.5e-06   &  \textbf{1.5e-06}   \\
31/01/2011   &  29.2   &  29.2  & 9.2e-06   &  \textbf{1.2e-06}   \\
10/05/2011   &  28.8   &  27.6  & 2.5e-06   &  \textbf{1.2e-06}   \\
10/06/2011   &  26.0   &  30.7  & \textbf{5.0e-07}   &  2.8e-06   \\
30/12/2011   &  20.3   &  30.0  & 6.8e-06   &  \textbf{3.6e-06}   \\
\hline
\end{tabular}
\end{table}

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GaussLOO.pdf}
%\subfloat[]{\label{GaussLOO}\includegraphics[scale=.4]{figures/GaussLOO.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/MatLOO.pdf}
%\subfloat[]{\label{MatLOO}\includegraphics[scale=.4]{figures/MatLOO.pdf}}
\end{minipage}
\caption{The function to be optimized in the LOO criterion (\ref{LOOLConst}) using the Gaussian (left) and the Mat\'ern 5/2 covariance function (right). Swap versus Euribor 6M on $30/12/2011$.}
%(resp. Mat\'ern 3/2) covariance function Figure~\ref{GaussLOO} (resp. Figure~\ref{MatLOO}).}
\label{LOOMatGauss}
\end{figure}

Once the length parameter $\theta$ has been estimated, the standard deviation parameter $\sigma$ is estimated using Equation~\eqref{sigSCV_2}. For instance, at quotation date $30/12/2011$, we obtain for respective Gaussian and Mat\'ern 5/2 covariance kernel $\hat{\sigma}_G=2.89$ and $\hat{\sigma}_M=0.93$.
%In Table~\ref{TABLETROIS}, we give estimated values of the variance parameter at each considered quotation date.
%\begin{table}[hptb]
%\centering
%\caption{\did{Variance parameters} estimation using ACV methods (Swap versus Euribor 6M).}
%\label{TABLETROIS}
%\begin{tabular}{ccc}
%\hline \hline
%Date  &  $\hat{\sigma}_G$ &  $\hat{\sigma}_M$ \\%& Gaussian optimal value & Mat\'ern 3/2 optimal value  \\   
%\hline
%02/06/2010   & 0.26 & 0.52 \\ %25.8   &  47.9  & 5.9e-05   &  \textbf{6.4e-05}   \\
%05/07/2010   &  & \\ %26.2   &  40.5  & 5.6e-05   &  1.1e-04   \\
%03/08/2010   &  & \\ %26.6   &  45.2  & 6.2e-05   &  8.8e-05   \\
%29/11/2010   &  & \\ %23.5   &  47.0  & 8.0e-06   &  8.9e-05   \\ 
%30/12/2010   &  & \\ %27.5   &  39.5  & 6.7e-05   &  9.1e-05   \\
%31/01/2011   &  & \\ %29.2   &  40.0  & 1.4e-04   &  8.5e-05   \\
%10/05/2011   &  & \\ %28.8   &  37.6  & \textbf{3.8e-05}   &  7.0e-05   \\
%10/06/2011   &  & \\ %20.0   &  36.5  & 5.2e-05   &  6.5e-05   \\
%30/12/2011   &  & \\ %20.3   &  33.0  & 1.0e-04   &  1.6e-04   \\
%\hline
%\end{tabular}
%\end{table}


%In Figure~\ref{LOOMatGauss}, we use for example the data correspond to Swap versus Euribor 6M on $02/06/2010$. Figure~\ref{GaussLOO} shows the function to be optimized in LOO criterion using the Gaussian covariance function, where Figure~\ref{MatLOO} shows the case using the Mat\'ern 3/2 covariance function. In Figure~\ref{GaussLOO}, the global optimum of the LOO criterion is complicated to find and the values of the function reaches to $0.7$ for some values of $\theta$ contrarily to Figure~\ref{MatLOO}, where the global optimal is simple and the values of the function are less than $0.05$ everywhere. \\


%\begin{figure}[hptb]
%\begin{minipage}{.5\linewidth}
%\centering
%\subfloat[]{\label{Gauss14}\includegraphics[scale=.4]{figures/Gauss14.pdf}}
%\end{minipage}%
%\begin{minipage}{.5\linewidth}
%\centering
%\subfloat[]{\label{Mat14}\includegraphics[scale=.4]{figures/Mat14.pdf}}
%\end{minipage}
%\caption{Simulated paths taken from the conditional GP with monotonicity (non-increasing) constraints and without the second to last information where $x=30$. The Gaussian (resp. Mat\'ern 3/2) kernel is used Figure~\ref{Gauss14} (res. Figure~\ref{Mat14}). The proposed estimator mode without the same information (black solid line). \hassan{Swap vs Euribor 6M data, quotation date 02/06/2010}.}
%\label{MatGauss14}
%\end{figure}
%
%
%In Figure~\ref{MatGauss14}, \hassan{we choose $N=40$ and generate 100 sample paths taken from model~\eqref{propModel} conditionally to given Euribor 6M data on $02/06/2010$ when removing the second to last \hassan{row} in matrix $A$ and vector $\boldsymbol{b}$ and monotonicity (non-increasing) constraints. It illustrates the increase of uncertainty when one observation is removed. This kind of prediction with removed observations is used in adapted cross validation techniques, for example to compute $\hat{c}^2_{i,-i}$ in Equation~\eqref{sigSCV}.} In Figure~\ref{Gauss14}, the Gaussian covariance function is used with \hassan{the parameters} $(\sigma,\hat{\theta}_G)=(1,25.8)$, where the Mat\'ern 3/2 one is used in Figure~\ref{Mat14} with \hassan{the associated parameters} $(\sigma,\hat{\theta}_M)=(1,47.9)$. In both figures, the black solid line represents the mode of the conditional GP without the second to last information $M_{i,\theta}(y_{-i})$, where $i=14$. Notice that it is non-increasing \hassan{in the entire domain} and respects the \hassan{linear equality condition (\ref{EqCond})}.

\paragraph{Curve construction at a single quotation date.}
We now illustrate the curve construction method described in Section \ref{ODC} in a one-dimensional setting. The construction is based on market quotes as of 
% the construction of term-structure of zero-coupon swaps rates
$30/12/2011$. For this particular date, the estimated length parameter for the Gaussian and the Mat\'ern 5/2 kernel are given in Table~\ref{EPOV} ($\hat{\theta}_G=20.3$ and $\hat{\theta}_M=30.0$). In that case, using the method described in Section~\ref{PE}, the estimated variance parameter is equal to $\hat{\sigma}_G=2.89$ for the Gaussian kernel and to $\hat{\sigma}_M=0.93$ for the Mat\'ern 5/2 kernel. 
 %Equation~(\ref{sigSCV}), we obtain $\hat{\sigma}_G=0.48$ which correspond to the variance parameter of the Gaussian covariance function and $\hat{\sigma}_M=0.93$ for the Mat\'ern 3/2 covariance function. \\
Figure~\ref{GaussMatParEst} compares  the sample paths of discount factors for the Gaussian and the Mat\'ern 5/2 covariance function using the corresponding estimated parameters. In both cases, we generate 100 sample paths taken from model~\eqref{propModel} conditionally  to  linear equality constraints~\eqref{EqXiCond} and non-increasing constraints. 
%Even if this is not visible in the graphic, the curves generated from the Gaussian covariance function are of class 
%%The Gaussian covariance function is used with the associated parameters (Figure~\ref{GaussParEst}) and the Mat\'ern 3/2 covariance function with the associated parameters (Figure~\ref{MatParEst})}. 
Note that the simulated curves (gray lines) are non-increasing in the entire domain. 
%\ares{Hassan: ce n'est pas visible sur le graphique
%As expected, these paths are smoother in Figure~\ref{GaussParEst} than  in Figure~\ref{MatParEst}. This is related to the smoothness of the Gaussian covariance function. 
%}
Additionally, the black solid line represents the most likely curve, i.e., the mode of the conditional GP. Recall that, by construction, this curve satisfies the given constraints. The black dashed-lines represent the 95\% point-wise confidence intervals quantified by simulation. 
Figure~\ref{SRGaussMatEURO} and Figure~\ref{TFGaussMatEURO} give the corresponding spot rates and instantaneous forward curves.\\

In order to compare our results with some models commonly used in most central banks, all figures are given together with the associated best-fitted Nelson-Siegel curves  \citep[see][]{nelson1987parsimonious} and the associated best-fitted Svensson curves \citep[see][]{Svensson:1994}. Parameters have been estimated by minimizing the sum of squared errors between market and model prices. We use a gradient descent algorithm with randomly chosen starting values as described in \cite{gilli2010calibrating}. The optimal parameters are given in Table~\ref{ParamNS_SWAP}. Discount factors and forward rates have been deduced from Nelson-Siegel and Svensson yield curves.\\

%Pour le jeu de données Swap vs Euribor 6M au 31/12/2011
%NS :
%lambda = 7.4615
%beta1 = 0.0189
%beta2 = -0.0160
%beta3 = 0.0487
%NSS :
%lambda1 = 4.0486
%lambda2 = 28.4285
%beta1 = 0.1719
%beta2 = -0.1590
%beta3 = -0.1101
%beta4 = -0.4093
\begin{table}[hptb]
\centering
\caption{Parameters estimation for Nelson-Siegel and Nelson-Siegel-Svensson model (Swap versus Euribor 6M on 30/12/2011).}
\label{ParamNS_SWAP}
\begin{tabular}{l|cccccc}
\hline \hline
 & $\lambda_1$ & $\lambda_2$ & $\beta_1$ & $\beta_2$ & $\beta_3$ & $\beta_4$ \\ 
\hline 
Nelson-Siegel          & 7.4615 &  - & 0.0189 &  -0.0160 & 0.0487 & - \\ \hline 
Nelson-Siegel-Svensson & 4.0486 & 28.4285 & 0.1719 & -0.1590 & -0.1101 & -0.4093 \\ 
\hline 
\end{tabular}
\end{table}




\begin{figure}[H]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.37]{figures/GaussParEst.pdf}
%\subfloat[]{\label{GaussParEst}\includegraphics[scale=.37]{figures/GaussParEst.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.37]{figures/DFNSswap.pdf}
\includegraphics[scale=.37]{figures/MatParEst.pdf}
%\subfloat[]{\label{MatParEst}\includegraphics[scale=.37]{figures/MatParEst.pdf}}
\end{minipage}
\caption{Simulated paths (gray lines) taken from the conditional GP with non-decreasing constraints and market-fit constraints using the Gaussian covariance function with nugget equal to $10^{-5}$ (left) and the Mat\'ern 5/2 covariance function without nugget (right).
%The simulated paths and the mode respect the linear equality condition defined in (\ref{EqXiCond}) and they are non-increasing \hassan{in the entire domain}. 
Swap vs Euribor 6M market quotes as of $30/12/2011$.}
\label{GaussMatParEst}
\end{figure}

\begin{figure}[H]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.37]{figures/SRGaussienEURO.pdf}
%\subfloat[]{\label{SRGaussienEURO}\includegraphics[scale=.37]{figures/SRGaussienEURO.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.37]{figures/SRNSswap.pdf}
\includegraphics[scale=.37]{figures/SRMatern32EURO.pdf}
%\subfloat[]{\label{SRMatern32EURO}\includegraphics[scale=.37]{figures/SRMatern32EURO.pdf}}
\end{minipage}
\caption{Spot rates obtained from sample paths of Figure \ref{GaussMatParEst} with Gaussian covariance function (left) and Mat\'ern 5/2 covariance function (right). Gray lines represent $-\frac{1}{x}\log Y^N(x)$ for each sample path. The black solid line is the most likely spot rate curve $-\frac{1}{x}\log M^N_K\left(x \suchthat A,\boldsymbol{b}\right)$.}
\label{SRGaussMatEURO}
\end{figure}

\begin{figure}[H]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.37]{figures/TFGaussienEURO.pdf}
%\subfloat[]{\label{TFGaussienEURO}\includegraphics[scale=.37]{figures/TFGaussienEURO.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.37]{figures/TFNSswap.pdf}
\includegraphics[scale=.37]{figures/TFMatern32EURO.pdf}
%\subfloat[]{\label{TFMatern32EURO}\includegraphics[scale=.37]{figures/TFMatern32EURO.pdf}}
\end{minipage}
\caption{Forward rates obtained from sample paths of Figure \ref{GaussMatParEst}  with Gaussian covariance function (left) and Mat\'ern 5/2 covariance function (right). Gray lines represent $-\frac{d}{dx}\log Y^N(x)$ for each sample path. The black solid line is the most likely forward rate curve $-\frac{d}{dx}\log M^N_K\left(x\suchthat A,\boldsymbol{b}\right)$.} 
%Gaussian (resp. Matérn 3/2) covariance function is used in Figure~\ref{TFGaussienEURO} (resp. Figure~\ref{TFMatern32EURO}). }
\label{TFGaussMatEURO}
\end{figure}

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[trim={0 5cm 0 0}, scale=.37]{figures/renteGaussSwap.pdf}
%\subfloat[]{\label{renteGaussSwap}\includegraphics[scale=.4]{figures/renteGaussSwap.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[trim={0 5cm 0 0}, scale=.37]{figures/renteMatSwap.pdf}
%\subfloat[]{\label{renteMatSwap}\includegraphics[scale=.4]{figures/renteMatSwap.pdf}}
\end{minipage}
\caption{ Histogram of a periodic annuity-due present values constructed from 100 000 independent simulations of discount factor curves under a Gaussian covariance function (left) and under a  Matérn covariance function (right). 
}
\label{renteGaussMatSwap}
\end{figure}

It can be seen on Figures~\ref{GaussMatParEst},~\ref{SRGaussMatEURO} and~\ref{TFGaussMatEURO} that the Nelson-Siegel model and the Nelson-Siegel-Svensson model are close to the monotonic kriging mode, this is especially the case for the Svensson extension. However, these models do not always fulfil the given market constraints (especially for the Nelson-Siegel model). Furthermore, they do not provide confidence intervals.

\begin{rem}
\label{rem:comparison_kernels}
As can be seen on Figure~\ref{TFGaussMatEURO}, the shape of the 
most probable forward curves (black solid lines) differs  under the two considered covariance kernels. When using a Gaussian covariance kernel, the mode curve features an extraneous hump for maturities greater than 30 years  whereas, in the Mat\'ern 5/2 case, it is slightly increasing and follows the fitted Nelson-Siegel-Svensson curve. This difference comes from the fact that $\mathcal{C}^{\infty}$ forward curves generated by the Gaussian kernel are  less flexible to adjust to local tendency of the data than  $\mathcal{C}^{1}$ forward curves generated by the Mat\'ern 5/2 kernel. Note moreover that, in the Gaussian kernel case, inverting the covariance matrix without a nugget effect  may involve numerical instabilities due to the restriction of infinite differentiability. This issue has also been reported in \cite{doi:10.1137/140976741}, where the authors state in particular that ``\textit{this choice of the Mat\'ern covariance function over the commonly used squared exponential family (Sacks et al., 1989) avoids numerical instability, often observed when inverting the covariance matrix, by removing the restriction of infinite differentiability''}.

 %This problem can be solved by introducing a small nugget effect.  
% This also explains the undesirable oscillations of confidence bands in the Gaussian kernel case. 
% This conclusion is in line with the empirical investigation of 
%\cite{doi:10.1137/140976741}. 
% Golchi et al. 2015 (SIAM/ASA JUQ) et dans cet article voir la phrase: 'This choice of the Matern covariance function over the commonly used squared exponential family (Sacks et al., 1989) avoids numerical instability' 
% Moreover, in the Mat\'ern 5/2  case,  the mode curve is closer to the pointwise mean of the sampled forward curves 
 \end{rem}


The previous simulations can be used to estimate the distribution of other  financial assets whose values depend on the curve. As an example, in Figure~\ref{renteGaussMatSwap}, we plot an histogram of the present value $\text{\"a}_{\annu{n}}^{(p)}=\sum_{k=0}^{pn-1}\frac{1}{p}Y^N(k/p)$ of a periodic annuity-due using 100 000 simulations of discount factors sample paths, where $n=40$ and $p=12$ (monthly payments). Notice that, despite the variability of the simulated sample paths of the conditional Gaussian process, the present value of the periodic annuity-due remains stable with 95\% confidence intervals $[25.12, 25.30]$ using the Gaussian kernel and $[25.07, 25.41]$ using the Mat\'ern kernel.%
% Figure~\ref{renteGaussSwap}
 %in Figure~\ref{renteMatSwap}.

% Gaussian
%     2.5%      50%    97.5% 
%    22.88232 23.01425 23.13055 

% Matérn 3/2
%    2.5%      50%    97.5% 
%  22.92277 23.03023 23.20211 


\paragraph{Several quotation dates.}


We now illustrate the building procedure in dimension two, when  data observed at different quotation dates are incorporated. We then construct a surface representing the evolution of discount factors as a function of time-to-maturities and quotation dates. To do this, we use the approach described in Section~\ref{TDC}. In Figure~\ref{Dim2Gauss}, the  surface represents the mode estimator of the conditional GP in dimension two. 
%based on the quotation dates defined in Table~\ref{EPOV}
The construction relies on the swap quotations at the 9 dates given in Table \ref{EPOV}. We choose $N_x=40$ and $N_t=20$ and we consider a two-dimensional Gaussian kernel written as
%\begin{equation*}
%K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x_1-x_1')^2}{2\theta_1}-\frac{(x_2-x_2')^2}{2\theta_2}\right),
%\end{equation*}
\begin{equation*}
K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x-x')^2}{2\theta_1}-\frac{(t-t')^2}{2\theta_2}\right),
\end{equation*}
where $\boldsymbol{x}=(x,t)$ and $\boldsymbol{x}'=(x',t')$. For each vector $\boldsymbol{x}=(x,t)$, the first component $x$ represents a time-to-maturity and the second component $t$ represents a quotation date. Without loss of generality, the distance $t'-t$ between two quotation dates has been expressed in percentage of the length between the two extreme dates of the sample. The parameters $\theta_1$ and $\theta_2$ are fixed respectively to $25$ and $0.5$. Notice that the constructed discount factor surface is non-increasing with respect to time-to-maturities.
%If we fix the second variable at any dates, the corresponding non-increasing curve represents the discount factor at this date.
%where the first variable represents the evaluation of the discount factors and the second one for the different quotation dates. 
%\ares{Hassan: je ne comprends pas cette phrase: Without loss of generality, the second variables is supposed to belong in the unit interval $[0,1]$.}
%\begin{figure}[hptb]

\begin{figure}[H]
\centering
\includegraphics[scale=0.35]{figures/Swap2D.pdf}
\caption{Swap-vs-Euribor discount factors as a function of  time-to-maturities and quotation dates.}
\label{Dim2Gauss}
\end{figure}

\subsection{OIS discount curves}\label{sec:OIS}

We now apply the kriging method to construct OIS discount curves. The aim is to construct, at different quotation dates $t$, discount curves $T\rightarrow P^D(t,T)$ based on market quotes of Overnight Indexed Swaps associated with different standard maturities.
We consider par rates of OIS at the 10 quotation dates given in Table~\ref{EPOVOIS}. 
For each of these quotation dates, the term-structure is built  from $14$ swap rates $S_1, \ldots, S_{14}$ associated with standard maturities in the set $E:=\{1,\ldots,10,15,20,30,40\}^\top$. 
For each standard maturity $T\in E$, the value $P(t,k)$ of the curve at time horizons $k=1, \ldots, T$ are linked through the linear relation~\eqref{eq:OIS}. Then, the curve is compatible with market quotes if the vector of discount factors $P^D(t, X) := (P^D(t, 1), \ldots, P^D(t,40))^{\top}$ satisfies a linear system of the form 
\begin{equation}
\label{eq:market_fit_OIS}
A_t \cdot P(t, X)=\boldsymbol{b}_t,
\end{equation}
where $A_t$ is  a $14\times 40$ real matrix and $\boldsymbol{b}_t= (1,\ldots,1)^\top\in \mathbb{R}^{14}$. In this case, we have $n=14$ observations which depends on $m=40$ points of the curve. Note that the market fit condition has exactly the same form as for the previous example in Subsection~\ref{CD}.
\\


\paragraph{Parameters estimation.}
In Table~\ref{EPOVOIS}, we estimate the length hyper parameter $\hat{\theta}_G$ and $\hat{\theta}_M$ associated respectively to Gaussian and Mat\'ern 5/2 covariance functions (see Table~\ref{kernel}). The optimal values in the last two columns of Table~\ref{EPOVOIS} correspond to the values of the LOO criterion defined in (\ref{LOO}) at the global optimum. Note that the estimated parameters  $\hat{\theta}_G$ and $\hat{\theta}_M$ remain stable across the considered quotation dates.  In addition, the obtained optimal values  are close  for the two covariance functions, even if they are slightly smaller for the Gaussian covariance function. Figure~\ref{GaussMatLOOOIS} represents the function to be optimized in criterion (\ref{LOO}) using the OIS data on $03/06/2010$. Given the shape of the functions, the estimation procedure has turned to be much more straightforward using a Mat\'ern 5/2 covariance function.

% (Figure~\ref{MatLOOOIS}), 
%
% to Figure~\ref{GaussLOOOIS}, where the Gaussian covariance function is used. \\


% Note also that the minimal value of the objective function in the LOO criterion \eqref{LOOLConst} is slightly smaller when using a Gaussian covariance function. However, as can be seen in Figure~\ref{LOOMatGauss}, the objective function (cross-validation error)  reach $0.7$ for the Gaussian covariance function whereas it never goes above $0.05$ for the Mat\'ern 3/2 covariance function. Additionally, the global minimum is easier to find in the Mat\'ern 3/2 case.


\begin{table}[hptb]
\centering
\caption{Parameters estimation using ACV methods (OIS data).}
\label{EPOVOIS}
\begin{tabular}{ccccc}
\hline \hline
Date  &  $\hat{\theta}_G$ &  $\hat{\theta}_M$ & Optimal value Gaussian & Optimal value Mat\'ern 5/2  \\   
\hline
03/06/2010 & 26.2 & 19.1 & \textbf{2.5e-05} & 9.5e-05 \\
04/10/2010 & 27.8 & 20.6 & \textbf{1.2e-05} & 7.4e-05 \\
31/12/2010 & 26.1 & 18.7 & \textbf{2.6e-05} & 8.4e-05 \\
04/03/2011 & 28.2 & 19.0 & \textbf{1.5e-05} & {4.7e-05} \\
15/06/2011 & 27.3 & 18.2 & \textbf{1.2e-05} & 6.9e-05 \\
10/10/2011 & 26.5 & 23.8 & \textbf{1.3e-05} & 5.3e-05  \\
14/11/2011 & 26.1 & 23.8 & \textbf{1.4e-05} & 7.4e-05 \\
15/12/2011 & 25.8 & 24.2 & \textbf{1.9e-05} & 7.6e-05 \\
\hline
\end{tabular}
\end{table}



\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GaussLOOOIS.pdf}
%\subfloat[]{\label{GaussLOOOIS}\includegraphics[scale=.4]{figures/GaussLOOOIS.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/MatLOOOIS.pdf}
%\subfloat[]{\label{MatLOOOIS}\includegraphics[scale=.4]{figures/MatLOOOIS.pdf}}
\end{minipage}
\caption{The function to be optimized in LOO criterion (\ref{LOOLConst}) using the Gaussian covariance function (left) and the Mat\'ern 5/2 covariance function (right). OIS data on 03/06/2010.
%(Figure~\ref{GaussLOOOIS}) and the Mat\'ern 5/2 covariance function (Figure~\ref{MatLOOOIS}).
}
\label{GaussMatLOOOIS}
\end{figure}

As previously, once the length parameter $\theta$ has been estimated, the standard deviation parameter $\sigma$ is estimated using Equation~\eqref{sigSCV_2}.

%\begin{figure}[hptb]
%\begin{minipage}{.5\linewidth}
%\centering
%\subfloat[]{\label{Gauss14OIS}\includegraphics[scale=.4]{figures/Gauss14OIS.pdf}}
%\end{minipage}%
%\begin{minipage}{.5\linewidth}
%\centering
%\subfloat[]{\label{Mat14OIS}\includegraphics[scale=.4]{figures/Mat14OIS.pdf}}
%\end{minipage}
%\caption{Simulation paths taken from the conditional GP with monotonicity (non-increasing) constraints and without the second to last information where $x=30$. The Gaussian (resp. Mat\'ern 5/2) kernel is used Figure~\ref{Gauss14OIS} (res. Figure~\ref{Mat14OIS}). The proposed estimator mode without the same information (black solid line).}
%\label{MatGauss14OIS}
%\end{figure}
%
%
%In Figure~\ref{MatGauss14OIS}, we choose $N=40$ and generate 100 sample paths taken from model (\ref{propModel}) conditionally to monotonicity (non-increasing) constraints and to linear equality conditions~\eqref{EqXiCond} without the second to last information given by the OIS data on 15/12/2011. In Figure~\ref{Gauss14OIS}, the Gaussian covariance function is used with the parameters \hassan{$(\sigma,\hat{\theta}_G)=(1,25.8)$}, where the Mat\'ern 5/2 one is used in Figure~\ref{Mat14OIS} with the parameters \hassan{$(\sigma,\hat{\theta}_M)=(1,60.7)$}. In both figures, the black solid line represents the mode of the monotone GP without the second to last information $M_{i,\theta}(y_{-i})$, where $i=14$. Notice that it is non-increasing in the entire domain and respects the linear equality condition defined in Equation~(\ref{EqCond}).

\paragraph{One single quotation date.}
In Figure~\ref{GaussMatParEstOIS}, we choose $N=50$ and generate 100 sample paths of discount factors constructed from model~\eqref{propModel} when using a Gaussian covariance function (left graph) and a Mat\'ern 5/2 covariance function (right graph). All the curves are non-increasing with respect to time-to-maturities. In addition, they all are perfectly compatible with  OIS data as of $03/06/2010$. The Gaussian process hyper-parameters have been estimated by the ACV method described in Section~\ref{PE}.
The estimated hyper-parameters are given by $(\hat{\theta}_G, \hat{\sigma}_G)=(26.2,\, 4.24)$ when using a Gaussian covariance function and by $(\hat{\theta}_M, \hat{\sigma}_M)=(19.1,\,0.24)$ when using a Mat\'ern 5/2 covariance function.
%These parameters have been estimated by the ACV method (Figure~\ref{GaussParEstOIS}), where the Mat\'ern 5/2 covariance function is used with the parameters $(\hat{\sigma}_M,\hat{\theta}_M)=(0.8,24)$ which are estimated by the ACV method (Figure~\ref{MatParEstOIS}).
% In both figures, the proposed estimator mode (black solid line) respects the linear equality condition defined in~(\ref{EqXiCond}) and it is non-increasing in the entire domain.
%Figure~\ref{GaussMatParEst} compares  the sample paths of discount factors for the Gaussian and the Mat\'ern 3/2 covariance function using the corresponding estimated parameters. In both cases, we generate 100 sample paths taken from model~\eqref{propModel} conditionally  to  linear equality constraints~\eqref{EqXiCond} and non-increasing constraints. 
%%Even if this is not visible in the graphic, the curves generated from the Gaussian covariance function are of class 
%%%The Gaussian covariance function is used with the associated parameters (Figure~\ref{GaussParEst}) and the Mat\'ern 3/2 covariance function with the associated parameters (Figure~\ref{MatParEst}). 
%Note that the simulated curves (gray lines) are non-increasing in the entire domain. 
%%\ares{Hassan: ce n'est pas visible sur le graphique
%%As expected, these paths are smoother in Figure~\ref{GaussParEst} than  in Figure~\ref{MatParEst}. This is related to the smoothness of the Gaussian covariance function. 
%%}
The black solid line represents the most likely curve, i.e., the mode of the conditional GP. Recall that, by construction, this curve satisfies the given constraints. The black dashed-lines represent the 95\% point-wise confidence intervals quantified by simulation. Figure~\ref{SRGaussMat} and Figure~\ref{TFGaussMat} gives the corresponding spot rate and forward curves.\\

Again, all figures~\ref{GaussMatLOOOIS},~\ref{SRGaussMat} and~\ref{TFGaussMat} are given together with the associated best-fitted Nelson-Siegel curves  \citep[see][]{nelson1987parsimonious} and the associated best-fitted Svensson curves \citep[see][]{Svensson:1994}.  Parameters have been estimated by minimizing the sum of squared errors between market and model prices. We use a gradient descent algorithm with randomly chosen starting values as described in \cite{gilli2010calibrating}. The optimal parameters are given in Table~\ref{ParamNS_OIS}.
%NS :
%lambda = 1.0890
%beta1 = 0.0341
%beta2 = -0.0171
%beta3 = -0.0601
%NSS :
%lambda1 = 1.0938
%lambda2 = 15.1891
%beta1 = 0.0502
%beta2 = -0.0164
%beta3 = -0.1074
%beta4 = -0.0494
\begin{table}[hptb]
\centering
\caption{Parameters estimation for Nelson-Siegel and Nelson-Siegel-Svensson model (OIS data, 03/06/2010).}
\label{ParamNS_OIS}
\begin{tabular}{l|cccccc}
\hline \hline
 & $\lambda_1$ & $\lambda_2$ & $\beta_1$ & $\beta_2$ & $\beta_3$ & $\beta_4$ \\ 
\hline 
Nelson-Siegel          & 1.0890 &  - & 0.0341 &  -0.0171 & -0.0601 & - \\ 
\hline 
Nelson-Siegel-Svensson & 1.0938 & 15.1891 & 0.0502 & -0.0164 & -0.1074 & -0.0494 \\ 
\hline 
\end{tabular}
\end{table}


\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GaussParEstOIS.pdf}
%\subfloat[]{\label{GaussParEstOIS}\includegraphics[scale=.4]{figures/GaussParEstOIS.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.4]{figures/DFNS.pdf}
\includegraphics[scale=.4]{figures/MatParEstOIS.pdf}
%\subfloat[]{\label{MatParEstOIS}\includegraphics[scale=.4]{figures/MatParEstOIS.pdf}}
\end{minipage}
\caption{OIS discount factor curves (gray lines) given as simulated paths of a conditional GP with non-increasing constraints using a Gaussian covariance function with nugget equal to $10^{-5}$ (left) and a Mat\'ern 5/2 covariance function without nugget (right). OIS data of 03/06/2010.
%All the curves are perfectly compatible with OIS par rates as of 
%with monotonicity constraints using the Gaussian covariance function with the parameters $(\hat{\sigma}_G,\hat{\theta}_G)=(1,25.8)$ Figure~\ref{GaussParEst} and the Mat\'ern 5/2 covariance function with the parameters $(\hat{\sigma}_M,\hat{\theta}_M)=(0.8,24)$ Figure~\ref{MatParEst}. Notice that in both figures, the simulated paths and the estimator mode
%are non-increasing in the entire domain and respect the linear equality condition defined in~(\ref{EqCond}).
}
\label{GaussMatParEstOIS}
\end{figure}

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/SRGaussienOIS.pdf}
%\subfloat[]{\label{SRGaussien}\includegraphics[scale=.4]{figures/SRGaussienOIS.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.4]{figures/SRNS.pdf}
\includegraphics[scale=.4]{figures/SRMatern52OIS.pdf}
%\subfloat[]{\label{SRMatern32}\includegraphics[scale=.4]{figures/SRMatern52OIS.pdf}}
\end{minipage}
\caption{
Spot rates obtained from sample paths of Figure \ref{GaussMatParEstOIS} with Gaussian covariance function (left) and Mat\'ern 5/2 covariance function (right). Gray lines represent $-\frac{1}{x}\log Y^N(x)$ for each sample path. The black solid line is the most likely spot rate curve $-\frac{1}{x}\log M^N_K\left(x \suchthat A,\boldsymbol{b}\right)$.
%Spot rates for sample paths and mode of the conditional Gaussian process. Gray lines represent $-\frac{1}{x}\log Y^N(x)$ given linear and monotonicity constraints for different sample paths. The black line corresponds to $-\frac{1}{x}\log M^N_K\left(x \suchthat A,\boldsymbol{b}\right)$. Gaussian (resp. Matérn 5/2) covariance function is used in Figure~\ref{SRGaussien} (resp. Figure~\ref{SRMatern32}). 
}
\label{SRGaussMat}
\end{figure}

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[trim={0 5cm 0 0}, scale=.4]{figures/TFGaussien.pdf}
%\subfloat[]{\label{TFGaussienOIS}\includegraphics[scale=.4]{figures/TFGaussien.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=.4]{figures/TFNS.pdf}
\includegraphics[trim={0 5cm 0 0}, scale=.4]{figures/TFMatern52.pdf}
%\subfloat[]{\label{TFMatern32OIS}\includegraphics[scale=.4]{figures/TFMatern52.pdf}}
\end{minipage}
\caption{
Forward rates obtained from sample paths of Figure  \ref{GaussMatParEstOIS}   with Gaussian covariance function (left) and Mat\'ern 5/2 covariance function (right). Gray lines represent $-\frac{d}{dx}\log Y^N(x)$ for each sample path. The black solid line is the most likely forward rate curve $-\frac{d}{dx}\log M^N_K\left(x\suchthat A,\boldsymbol{b}\right)$.
%Forward rates for sample paths and mode of the conditional Gaussian process. Gray lines represent $-\frac{d}{dx}\log Y^N(x)$ given linear and monotonicity constraints for different sample paths. The black line corresponds to $-\frac{d}{dx}\log M^N_K\left(x \suchthat A,\boldsymbol{b}\right)$. Gaussian (resp. Matérn 5/2) covariance function is used in Figure~\ref{TFGaussienOIS} (resp. Figure~\ref{TFMatern32OIS}). 
}
\label{TFGaussMat}
\end{figure}

Regarding comparison of curves constructed by the two considered covariance kernels, Remark \ref{rem:comparison_kernels} also applies here.

It appears that Nelson-Siegel and Nelson-Siegel-Svensson models do not fulfil all market constraints, contrary to the proposed methodology, which also gives confidence intervals.

%We now illustrate the building procedure in dimension two, when  data observed at different quotation dates are incorporated. We then construct a surface representing the evolution of discount factors as a function of time-to-maturities and quotation dates. To do this, we use the approach described in Section~\ref{TDC}. In Figure~\ref{Dim2Gauss}, the  surface represents the mode estimator of the conditional GP in dimension two. 
%%based on the quotation dates defined in Table~\ref{EPOV}
%The construction relies on the swap quotations at the 9 dates given in Table \ref{EPOV}. We choose $N_x=40$ and $N_t=20$ and we consider a two-dimensional Gaussian kernel written as
%\begin{equation*}
%K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x-x')^2}{2\theta_1}-\frac{(t-t')^2}{2\theta_2}\right),
%\end{equation*}
%where the variable $x$ represents time-to-maturities and the variable $t$ represents  quotation dates. \ares{Without loss of generality,  the series of quotation dates has been translated into a series of values in the interval $[0,1]$, so that the distance unit between quotation dates $t-t'$ is expressed as a percentage in $[0,1]$ of the considered quotation period length. The parameters $\theta_1$ and $\theta_2$ are fixed respectively to $25$ and $0.5$. Notice that the constructed discount factor surface is non-increasing with respect to time-to-maturities.}

\paragraph{Several quotation dates.}
%The study at one single quotation date was relying on the one dimensional monotone interpolation technique described in Section~\ref{ODC}. 
Using the two dimensional approach described in Section~\ref{TDC}, we build in Figure~\ref{Dim2GaussOIS}  a surface representing OIS discount factors with respect to time-to-maturities and quotation dates. The construction relies on OIS quotations at the 8 dates given in Table \ref{EPOVOIS}. The surface corresponds to the mode of the conditional GP given market-fit equality constraints and non-increasing constraints in the direction of time-to-maturities.  We choose $N_x=40$ and $N_t=20$ and we consider a two-dimensional Gaussian kernel written as
%\begin{equation*}
%K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x-x')^2}{2\theta_1}-\frac{(t-t')^2}{2\theta_2}\right),
%\end{equation*}
%where the variable $x$ represents time-to-maturities and the variable $t$ represents  quotation dates.
\begin{equation*}
K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x-x')^2}{2\theta_1}-\frac{(t-t')^2}{2\theta_2}\right),
\end{equation*}
%where $\theta_1$ and $\theta_2$ are fixed respectively to $25$ and $0.5$.
where $\boldsymbol{x}=(x,t)$ and $\boldsymbol{x}'=(x',t')$. For each vector $\boldsymbol{x}=(x,t)$, the first component $x$ represents a time-to-maturity and the second component $t$ represents a quotation date.
Without loss of generality, the distance  $t-t'$ between two quotation dates has been expressed in percentage of the length between the two extreme dates of the sample. The parameters $\theta_1$ and $\theta_2$ are fixed respectively to $25$ and $0.5$. Observe that the constructed discount factor surface is non-increasing with respect to time-to-maturities.


%\begin{figure}[hptb]
\begin{figure}[H]
\centering
\includegraphics[scale=0.35]{figures/Dim2GaussOIS.pdf}
\caption{
OIS discount factors as a function of  time-to-maturities and quotation dates.
%The simulated surface represents the mode of the conditional GP with monotonicity constraints with respect to the first variable only using the OIS data which are defined in Table~\ref{EPOVOIS}.
}
\label{Dim2GaussOIS}
\end{figure}






\subsection{CDS implied default distributions}

We now apply the kriging method to build CDS implied default distributions. Indeed, CDS can be seen as an insurance product that covers the  loss of a particular debt issuer if the latter defaults over a certain protection period or maturity. The market quotes (CDS spreads) of these products provide information on the current cost for protection and, in turn,  on how the market assesses default probabilities of the underlying entity at different time horizons. 
\\

Our aim is to construct, at different quotation dates $t$, implied survival functions $T\rightarrow Q(t,T)$ of a particular debt issuer by observing the corresponding term-structure of CDS spreads (CDS spreads at increasing protection horizons). The quantity $Q(t,T)$ gives the probability at time $t$ that the underlying debt issuer does not default before the time horizon $T$. In this numerical illustration, we consider CDS  on the Russian sovereign debt for the 10 quotation dates given in Table~\ref{EPOVCDS}. 
For each of these quotation dates, this corresponds to $7$ CDS spreads $S_1, \ldots, S_{7}$ associated with protection maturities (in years) in the set $E:=\{1,2,3,4,5,7,10\}^\top$. 
For each standard maturity $T\in E$, the values of the survival probabilities at each premium payment dates $\ta_1<\cdots<\ta_n=T$ are linked through the linear relation~\eqref{eq:discrete_linear_CDS}. Note that the  premium payment dates are separated by a quarter period and these dates coincide for all quoted CDS maturities. 
In our numerical illustrations, the expected recovery rate $R$ is fixed to $40\%$ and the discount factors $P^D(t, \tau_k)$ are constructed by linear interpolation of treasury constant maturity rates given (for all considered quotation dates) by the Federal Reserve Bank of St. Louis. Then, if $\delta = 1/4$ represents a quarter period, the implied default distribution is compatible with market quotes if the vector of survival probabilities  $Q(t, X) := (Q(t, \delta), Q(t, 2\delta), \ldots, Q(t,10))^{\top}$ satisfies a linear system of the form 
\begin{equation}
\label{eq:market_fit_CDS}
A_t \cdot Q(t, X)=\boldsymbol{b}_t
\end{equation}
where $A_t$ is  a $7\times 40$ real matrix and $\boldsymbol{b}_t= (1-R,\ldots,1-R)^\top\in \mathbb{R}^{7}$. In this case, we have $n=7$ observations which depends on $m=40$ points of the curve.


\paragraph{Parameters estimation.}
In Table~\ref{EPOVCDS}, we compare the estimation of  the length hyper parameter  for a Gaussian kernel ($\hat{\theta}_G$) and for a Mat\'ern 5/2 kernel  ($\hat{\theta}_M$). The estimation has been performed using the ACV method described in Section~\ref{PE}. 
%Note that the  estimated parameters $\hat{\theta}_G$ (Gaussian kernel) and $\hat{\theta}_M$ (Mat\'ern 3/2 kernel) remain stable across all considered quotation dates. 
As can be seen on Figure \ref{LOOMatGaussCDS}, the LOO objective function looks similar for these two covariance functions. 
%Using CDS quotes as of $23/08/2012$, the optimal value of the LOO criterion is equal to $6.5\times 10^{-5}$ at $\hat{\theta}_G=4.0$ (Figure~\ref{GaussLOOCDS}) and it is equal to $4.3\times 10^{-6}$ at $\hat{\theta}_M=15.8$ (Figure~\ref{MatLOOCDS}). 

\begin{table}[hptb]
\centering
\caption{Length parameters estimation using ACV methods (CDS data).}
\label{EPOVCDS}
\begin{tabular}{ccccc}
\hline \hline
Date  &  $\hat{\theta}_G$ &  $\hat{\theta}_M$ & Optimal value Gaussian & Optimal value Mat\'ern 5/2   \\   
\hline
06/01/2005 & 6.6  & 10.5  & 6.3e-06 & \textbf{4.9e-07} \\
02/02/2006 & 4.9  & 15.5  & {2.7e-07} & \textbf{1.8e-07} \\
20/03/2007 & 4.5  & 19.8  & 7.6e-06 & \textbf{7.1e-07} \\
04/04/2008 & 9.7  & 11.4  & 5.0e-06 & \textbf{2.6e-06} \\
11/05/2009 & 4.9 & 10.1  & 1.3e-04  & \textbf{4.2e-05} \\
21/06/2010 & 11.1 & 10.1  & 1.1-05  & \textbf{3.4e-06}  \\
14/07/2011 & 14.5 & 17.1  & \textbf{1.2e-06} & 1.8e-06 \\
23/08/2012 & 4.1  & 10.8  & 8.2e-06 & \textbf{1.5e-06} \\
\hline
\end{tabular}
\end{table}

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/GaussLOOCDS.pdf}
%\subfloat[]{\label{GaussLOOCDS}\includegraphics[scale=.4]{figures/GaussLOOCDS.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=.4]{figures/MatLOOCDS.pdf}
%\subfloat[]{\label{MatLOOCDS}\includegraphics[scale=.4]{figures/MatLOOCDS.pdf}}
\end{minipage}
\caption{The function to be optimized in the LOO criterion (\ref{LOOLConst}) using the Gaussian (left) and the Mat\'ern 5/2 (right) covariance function. CDS quotes as of $06/01/2005$. 
%The optimal value of the LOO criterion is equal to $6.5\times 10^{-5}$ at $\hat{\theta}_G=4.0$ (Figure~\ref{GaussLOOCDS}) and it is equal to $4.3\times 10^{-6}$ at $\hat{\theta}_M=15.8$ (Figure~\ref{MatLOOCDS}). % Figure~\ref{GaussLOOCDS} (resp. Figure~\ref{MatLOOCDS}).
}
\label{LOOMatGaussCDS}
\end{figure}

%In Figure~\ref{LOOMatGaussCDS}, we use for example the data correspond to CDS on $23/08/2012$. Figure~\ref{GaussLOOCDS} shows the function to be optimized in LOO criterion~\eqref{LOOLConst} using the Gaussian covariance function. In Figure~\ref{MatLOOCDS} the Mat\'ern 3/2 covariance function is used. The optimal value of the LOO criterion is equal to $6.5\times 10^{-5}$ at $\hat{\theta}_G=4.0$ (Figure~\ref{GaussLOOCDS}) and it is equal to $4.3\times 10^{-6}$ at $\hat{\theta}_M=15.8$ (Figure~\ref{MatLOOCDS}). 

As previously, once the length parameter $\theta$ has been estimated, the standard deviation parameter $\sigma$ is estimated using Equation~\eqref{sigSCV_2}.

\paragraph{One single quotation date.}


In Figure~\ref{GaussMatCDS}, we choose $N=50$ and generate 100 sample paths of CDS implied survival curves, constructed from model~\eqref{propModel} when using a Gaussian covariance function (left graph) and a Mat\'ern 5/2 covariance function (right graph). All the curves are non-increasing with respect to time horizons. In addition, they all are perfectly compatible with  CDS data as of $06/01/2005$. 
The Gaussian process hyper-parameters have been estimated by the ACV method described in Section \ref{PE}.
The estimated hyper-parameters are given by  $(\hat{\theta}_G, \hat{\sigma}_G)=(4.9,\, 0.09)$ when using a Gaussian covariance function and by $(\hat{\theta}_M,\hat{\sigma}_M)=(10.5,\,0.22)$ when using a Mat\'ern 5/2 covariance function.
The black solid line represents the most likely curve, i.e., the mode of the conditional GP. Recall that, by construction, this curve satisfies the given constraints. The black dashed-lines represent the 95\% point-wise confidence intervals quantified by simulation. 

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[trim=0 5cm 0 0 0, scale=.4]{figures/GaussCDS.pdf}
%\includegraphics[scale=.4]{figures/GaussLOOCDS.pdf}
%\subfloat[]{\label{GaussCDS}\includegraphics[scale=.4]{figures/GaussCDS.pdf}}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[trim=0 5cm 0 0 0, scale=.4]{figures/MatCDS.pdf}
%\subfloat[]{\label{MatCDS}\includegraphics[scale=.4]{figures/MatCDS.pdf}}
\end{minipage}
\caption{
CDS implied survival curves (gray lines) given as simulated paths of a conditional GP with non-increasing constraints using a Gaussian covariance function (left)  or a  Mat\'ern 5/2 covariance function (right). 
%Simulated paths (gray lines) drawn from the conditional GP with monotonicity (non-increasing) constraints
%using the Gaussian (resp. Mat\'ern 3/2) covariance function Figure~\ref{GaussCDS} (resp. Figure~\ref{MatCDS}).
%Notice that the gray lines and the mode (black solid line) respect both equality conditions and monotonicity (non-increasing) constraints in the entire domain.
}
\label{GaussMatCDS}
\end{figure}

%As expected, in Figure~\ref{GaussCDS}, the simulated paths (gray lines) are more smoothing than ones in Figure~\ref{MatCDS}. This is because the smoothness of the Gaussian covariance function. In both figures, the proposed estimator mode (black solid line) which is defined in Equation~(\ref{propEstimator}) respects the linear equality condition defined in~(\ref{EqXiCond}) and the monotonicity (non-increasing) constraints in the entire domain. The black dashed-lines represent 95\% confidence intervals which are quantified from simulation.



\paragraph{Several quotation dates.}
Using the two dimensional approach described in Section~\ref{TDC}, we build in Figure~\ref{Dim2GaussCDS} a surface representing CDS implied survival curves as a function of time horizons and quotation dates.  The construction relies on CDS quoted spreads at the 8 dates given in Table~\ref{EPOVCDS}. 
The surface corresponds to the mode of the conditional GP given market-fit equality constraints and non-increasing constraints in the direction of time-to-maturities.  We choose $N_x=40$ and $N_t=20$ and we consider a two-dimensional Gaussian kernel written as
\begin{equation*}
K(\boldsymbol{x},\boldsymbol{x}')=\exp\left(-\frac{(x-x')^2}{2\theta_1}-\frac{(t-t')^2}{2\theta_2}\right),
\end{equation*}
%\did{In these numerical illustrations, the distance unit between quotation dates $t-t'$ is expressed as a percentage in $[0,1]$ of the considered quotation period length.}
where $\boldsymbol{x}=(x,t)$ and $\boldsymbol{x}'=(x',t')$. For each vector $\boldsymbol{x}=(x,t)$, the first component $x$ represents a time-to-maturity and the second component $t$ represents a quotation date.
Without loss of generality, the distance  $t-t'$ between two quotation dates has been expressed in percentage of the length between the two extreme dates of the sample. The parameters $\theta_1$ and $\theta_2$ are fixed respectively to $8$ and $1.7$. Notice that the constructed discount factor surface is non-increasing with respect to time-to-maturities. Considering several quotation dates simultaneously offers the advantage of increasing the data set for a better estimation of hyper-parameters  and of creating a consistent interpolation procedure across two directions (time horizons and quotation dates).

%where $\theta_1$ and $\theta_2$ are fixed respectively to $8$ and $1.7$. One can check that, in the entire domain, the resulting curve is decreasing with respect to the time-to-maturity only. Considering several quotation dates simultaneously offers the advantage of increasing the number of observations for hyper-parameters, to create a consistent interpolation for each quotation date, and to eventually consider non (yet) observed quotation dates.\\

%\begin{figure}[hptb]
\begin{figure}[H]
\centering
\includegraphics[scale=0.35]{figures/Dim2GaussCDS.pdf}
\caption{
CDS implied survival probabilities as a function of  time-to-maturities and quotation dates.
%The mode of the conditional GP with monotonicity constraints with respect to the first variable only using the CDS data defined in Table~\ref{EPOVCDS}.
}
\label{Dim2GaussCDS}
\end{figure}


\subsection{Need for monotonic kriging techniques}\label{sec:motivMonotonie}

%We try to interpolate the values of a studied quantity $Y(x)$ at any time $x$, given some market constraints. As seen before, a first set of observed market conditions can be summarized by
%\begin{equation*}
%A \mathbf{Y}(\mathbf{X})=\mathbf{b}\, ,
%\end{equation*}
%where $A$ is not necessarily a square matrix, so that these constraints do not lead to a unique set of possible values for $Y(x^{(i)})$, $i=1, \ldots, n$.

At the beginning of the Section~\ref{IntModIConst}, we justified monotonic techniques developed in this paper by some no-arbitrage constraints.
A natural question is whether a monotonic condition is required to construct realistic term-structures. Let us stress that the construction problems we have considered rely on a market information summarized by an ill-posed system of equations of the form
\begin{equation*}
A \mathbf{Y}(\mathbf{X})=\mathbf{b}\, ,
\end{equation*}
where $A$ is not necessarily a square matrix. Obviously, these constraints do not lead to a unique set of possible values for $Y(x^{(i)})$, $i=1, \ldots, m$, so that spot rates or implied default probabilities are not directly available at point $x^{(1)}, \ldots, x^{(m)}$.
In the absence of arbitrage opportunities, the process $Y$ representing default-free zero-coupon bonds, discount factors or implied survival probabilities shall be monotone. It is thus natural to consider a constrained interpolation technique to construct curves based on such quantities.\\

Consider here the previous swap curve construction problem but relax the monotonicity constraint on the default-free zero-coupon bond process $Y$. As can be seen on Figure~\ref{fig:motivMonotoneDF} (left panel), the resulting \emph{kriging mean  is not  a decreasing function}. In addition, even when the average curve is monotonic, some \textit{sample curves are clearly not monotonic}, leading to wide and \textit{unrealistic confidence intervals}. On the right panel of Figure~\ref{fig:motivMonotoneDF} one can see that corresponding spot rates may seem reasonable in average, but omitting the monotonicity constraints leads to wide confidence intervals.
%\emph{confidence intervals are too wide} or unsuited to monotonicity constraints.
%We motivate here the need of specific monotone interpolation techniques instead classical non constraint ones, thus justifying techniques that are developed in further Sections~\ref{ODC} and~\ref{TDC}. 
%We first discuss the situation where monotonic quantities are chosen to be directly interpolated. Then, we consider  the situation where non-monotone quantities such as spot rates are chosen to be directly interpolated.
%\paragraph{Interpolation of monotone quantities}\label{sec::motiv_monotoneCase}
%When trying to interpolate directly the monotone curve $Y$, it is natural to use monotonic interpolation techniques. 
%Two augments go in favor of direct interpolation of monotonic quantities.
%Non monotonic interpolation techniques are sometimes used in the Literature \did{A COMPLETER. Citer les références du REVIEWER 3 et 2}.%Here are the references mentioned by the reviewers.
%\begin{itemize}
%\item Reviewer 2 : \cite{LauriniOhashi2015}, \cite{Fengler:2015}, \cite{Laurini:2010}
%\item Reviewer 3 : \cite{Kanevski2008}, \cite{Benth:2015}, \cite{Steeley:2008}, \cite{nelson1987parsimonious}, \cite{Svensson:1994}
%\end{itemize}
%In this setting, it seems necessary to use monotonic interpolation techniques:
%\begin{itemize}
%\item \emph{No-Arbitrage constraints} lead to a theoretical need of monotonicity in curve constructions. The relationship of shape and monotonicity restrictions with no-arbitrage conditions is detailed for example in the recent literature, see~\cite{Fengler:2015} and~\cite{Laurini:2010}. Our setting thus corresponds to \emph{classical constrained setting}, which appears e.g. in constrained splines literature.\did{ARESKI: en dire un peu plus?}
%\item Without monotonicity constraints, the resulting \emph{average interpolation is not necessarily monotonic}. Furthermore, even when the average curve is monotonic, \emph{confidence intervals are too wide} or unsuited to monotonicity constraints.
%\end{itemize}
%\didsuppr{These problems are illustrated in Figure~\ref{fig:motivMonotoneDF}. On the left panel of this figure, one can see a case where the kriging mean is locally increasing and thus unsuited to discount factors. More interestingly, on the right panel, one can see that even when the kriging mean is decreasing, confidence intervals are very wide and correspond to non-monotonic paths, thus unsuited to discount factors.}
Notice that Figure~\ref{fig:motivMonotoneDF} does just aim at illustrating some problems that can be encountered when using unconstrained techniques, so that other choices of kernel or parameters in an unconstrained setting are not discussed here.\\

\begin{figure}[hptb]
\begin{minipage}{.5\linewidth}
\centering
%\includegraphics[scale=0.4]{figures/motivMonotoneDF2.pdf}
\includegraphics[scale=0.4]{figures/unconstrainedModel1e-7.pdf}
\end{minipage}%
\begin{minipage}{.5\linewidth}
\centering
\includegraphics[scale=0.4]{figures/unconstrainedModel1e-5.pdf}
\end{minipage}
%\caption{\did{Interpolation of some discount factors obtained with standard kriging techniques without monotonicity constraints, under a linear constraint $A \mathbf{Y}(\mathbf{X})=\mathbf{b}$ (the difference between the two curves relies on the choice of different length correlation parameters for the covariance kernel, the data is the one described at section~\ref{sec:OIS}).}}
\caption{Interpolation of some discount factors (left panel) using standard kriging techniques without monotonicity constraints. OIS Data, 03/06/2010, Gaussian kernel with nugget $10^{-7}$%(left panel) or $10^{-5}$ (right panel).
. Corresponding spot rates are given in the right panel.
In this figure only, linear constraints $A \mathbf{Y}(\mathbf{X})=\mathbf{b}$ has been approximated by point-wise constraints $Y(x^{(i)})= M_K^N(x^{(i)}|A,\mathbf{b})$, $i=1, \ldots, m$, in order to use standard unconstrained R package \textit{DiceKriging} \citep[see][]{Roustant:Ginsbourger:Deville:2012:JSSOBK:v51i01}. Estimated parameters from this package are $(\theta,\sigma^2)=(5.8, 0.22)$.}
\label{fig:motivMonotoneDF}
\end{figure}

%\paragraph{Direct interpolation of spot rates or forward rates}\label{sec::motiv_nonMonotoneCase}

If one wants to avoid using monotone interpolation techniques, a natural idea is to interpolate some quantities that are deduced from $Y(x)$ and that are not necessarily monotonic.
As an example, one can work on a function $\zeta(x)$ deduced in a bijective way from each path of $Y(x)$. For example in a financial setting, if $Y(x)$ are discount factors, the function $\zeta(x)$ could be spot rates, $\zeta(x)= -\frac{1}{x}\log Y(x)$, or forward rates, $\zeta(x)= -\frac{d}{dx} \log Y(x)$. Obviously these functions does not suffer from monotonicity constraints, so that this could avoid using monotone interpolation. In the literature, some studies interpolate interest rates without considering any monotonic constraints. \cite{Steeley:2008} concludes that it is better to directly interpolates spot rates rather than discount factors: ``\textit{better yield curves estimates are obtained by fitting to the yield curve directly rather than fitting first to the discount function}''. %However, this conclusion is based on observed market data consisting of sport 
 Some authors also propose to use kriging techniques to some non-monotonic curves, as~\cite{Benth:2015} who applies kriging to futures curves from energy futures prices. In~\cite{Kanevski2008}, tools from spatial statistics and machine learning are applied to produce some interest rate mapping  in
a two dimensional feature space (maturity, time). However, in these studies, rates or prices to be interpolated are assumed to be directly observed. In our setting, we do not necessarily observe  spot rates or even discount factors, and one aims at fitting market data in a non-parametric setting. We can make general objections to the interpolation of non-monotonic deduced quantities : %in the presented setting:
\begin{itemize}
\item The first objection is that even if $\zeta(x)$ is not monotone, it is \emph{still constrained}: as an example, quantities such as discount rates or forward rates  are expected to be positive in absence of arbitrage opportunity. As for the discount rates, interpolating spot rates could lead to locally negative average spot rates, and confidence intervals obtained by kriging can reach the threshold 0, which is not desirable: the constraints are translated from monotonicity constraints to positivity constraints, and classical kriging techniques cannot handle these constraints.
\item A second objection is that in our setting, observations are $A \mathbf{Y}(\mathbf{X})=\mathbf{b}$, where $A$ is not necessarily a square and invertible matrix. Thus \emph{we cannot directly observe}, with our data, spot rates or forward rates. Even in the very special case where $A$ is a square and invertible matrix, spot rates can be deduced from $\mathbf{Y}(\mathbf{X})$, but not forward rates: knowing a function at some abscissas does not give straightforward constraints on its derivative.
\item At last, using kriging interpolation, even if omitting the positivity of $\zeta(x)$, would lead to \emph{non linear conditions} in terms of $\zeta(X)$. One has seen that conditionally to $A \mathbf{Y}(\mathbf{X})=\mathbf{b}$, $Y(x)$ is still a Gaussian Process, but it would not be the case any more for the process $\zeta(x)$.
Even in the simple case where $A$ is invertible, the process $\zeta(x)$ must be positive and given $A \mathbf{Y}(\mathbf{X})=\mathbf{b}$, it is not Gaussian any more, so that suited interpolation techniques still have to be introduced.
\end{itemize}

%In the Figure~\ref{fig:motivMonotoneSR}, we have deduced spot rates $\zeta(x)$ from paths of the discount rates of the right panel of Figure~\ref{fig:motivMonotoneDF}. One can see that as discount factors paths are not monotonic, they can reach low values (or even negative values) that induce very high unrealistic spot rates values (or even undefined values). When discount rates become non-positive, this naturally cause numerical problems (in this illustration treated with an absolute value over a small minimal value).
%\begin{figure}[hptb]
%\begin{minipage}{.9\linewidth}
%\centering
%\includegraphics[scale=0.3]{figures/motivMonotoneSR1.pdf}
%\end{minipage}%
%%\begin{minipage}{.5\linewidth}
%%\centering
%%\includegraphics[width=7cm, height=7cm]{figures/motivMonotoneSR2.pdf}
%%\end{minipage}
%\caption{Interpolation of some spot rates obtained with standard kriging techniques without positivity constraints, deduced from paths of unconstrained discount factors (right).}
%\label{fig:motivMonotoneSR}
%\end{figure}


\section*{Conclusion}
In this paper, we show how suitable kriging techniques can be used to quantify model uncertainty embedded in the construction of financial term-structures. We consider that the curve under construction is an unobservable path of a conditional (spatial) Gaussian process satisfying some linear equality constraints and monotone properties. A suitable cross-validation method is proposed to estimate the Gaussian process covariance  parameters that control the level of uncertainty. We then investigate the efficiency of the proposed approach on some illustrative examples in one and two dimensions. The generated curves are all compatible with market quotes and respect non-arbitrage conditions. The conditional Gaussian process also allows to derive confidence bands for financial term-structures and related quantities. We compare the Gaussian  and the Mat\'ern 5/2 covariance kernels over different data sets : swaps vs Euribor, overnight indexed swaps, credit default swaps. We conclude that, for these applications, the Mat\'ern 5/2 covariance kernel seems to be the more appropriate since it generates more realistic forward curves in comparison with fitted  Nelson-Siegel or Svensson models.
In this work, we do not fully investigate the impact of curve uncertainty on the assessment of related  products and their hedging strategies. Moreover, in the illustrative example we have presented, we  consider that  market information is observed without uncertainty. It may be the case  that, due to  lack of liquidity, market quotes cannot be considered to be reliable. The kriging techniques we have proposed could then be adapted to the presence of noisy observations  (see section \ref{subsec:noisy}).
 %using the so-called \textit{nugget effect}. 
Improvements on the rejection sampling algorithm could also be useful, especially for the estimation of hyper-parameters in dimension $2$ where a large number of quotation dates has to be considered. These points are leaved for future research.

\section*{Acknowledgement}
We would like to thank the three anonymous reviewers and the editor for the time spent on this article and for useful suggestions. The authors also thank Xavier Bay (EMSE) and Nicolas Durrande (EMSE) for helpful discussions. This work has been conducted within the frame of the ReDice Consortium,
gathering industrial (CEA, EDF, IFPEN, IRSN, Renault) and academic (Ecole des Mines de Saint-Etienne, INRIA, and the University of Bern) partners around advanced methods for Computer Experiments. It also benefits from the support of the GRI in Financial Services and the Louis Bachelier Institute.  One author also thanks the ANR Lolita research project and DAMI research project.
\bibliography{bibliographie}
\bibliographystyle{apalike}

\end{document}
\end{document}

Kriging_of_Financial_Term_Structure_06_avril_2016.bbl
\begin{thebibliography}{}

\bibitem[Abrahamsen, 1997]{abrahamsen1997review}
Abrahamsen, P. (1997).
\newblock {\em {A review of Gaussian random fields and correlation functions}}.
\newblock Norsk Regnesentral/Norwegian Computing Center.

\bibitem[Abrahamsen and Benth, 2001]{abrahamsen2001kriging}
Abrahamsen, P. and Benth, F.~E. (2001).
\newblock Kriging with inequality constraints.
\newblock {\em Mathematical Geology}, 33(6):719--744.

\bibitem[Ametrano and Bianchetti, 2009]{Ametrano2009}
Ametrano, F. and Bianchetti, M. (2009).
\newblock {\em {Bootstrapping the illiquidity: Multiple yield curves
  construction for market coherent forward rates estimation}}, chapter~1.
\newblock Risk Books.

\bibitem[Andersen, 2007]{Andersen2007}
Andersen, L. (2007).
\newblock Discount curve construction with tension splines.
\newblock {\em Review of Derivatives Research}, 10(3):227--267.

\bibitem[Asgharian et~al., 2013]{Asgharian20134738}
Asgharian, H., Hess, W., and Liu, L. (2013).
\newblock A spatial analysis of international stock market linkages.
\newblock {\em Journal of Banking and Finance}, 37(12):4738 -- 4754.

\bibitem[Bachoc, 2013]{Bachoc201355}
Bachoc, F. (2013).
\newblock {Cross Validation and Maximum Likelihood estimations of
  hyper-parameters of Gaussian processes with model misspecification}.
\newblock {\em Computational Statistics \& Data Analysis}, 66(0):55 -- 69.

\bibitem[Barzanti and Corradi, 1999]{Barzanti-Corradi99}
Barzanti, L. and Corradi, C. (1999).
\newblock A note on direct term structure estimation using monotonic splines.
\newblock {\em Rivista di matematica per le scienze economiche e sociali},
  22:101--108.

\bibitem[Bay et~al., 2015]{bay:hal-01136466}
Bay, X., Grammont, L., and Maatouk, H. (2015).
\newblock {A New Method for Interpolating in a Convex Subset of a Hilbert
  Space}.
\newblock hal-01136466.

\bibitem[{Bay} et~al., 2016]{2016arXiv160202714B}
{Bay}, X., {Grammont}, L., and {Maatouk}, H. (2016).
\newblock {Generalization of the Kimeldorf-Wahba correspondence for constrained
  interpolation}.
\newblock {\em ArXiv e-prints}.

\bibitem[Baysal et~al., 2008]{baysal2008response}
Baysal, R.~E., Nelson, B.~L., and Staum, J. (2008).
\newblock Response surface methodology for simulating hedging and trading
  strategies.
\newblock In {\em Simulation Conference, 2008. WSC 2008. Winter}, pages
  629--637. IEEE.

\bibitem[Benth, 2015]{Benth:2015}
Benth, F.~E. (2015).
\newblock Kriging smooth futures curves.
\newblock {\em Energy Risk}, February:64--69.

\bibitem[Branger and Schlag, 2004]{Branger2004}
Branger, N. and Schlag, C. (2004).
\newblock {Model risk: A conceptual framework for risk measurement and
  hedging}.
\newblock In {\em EFMA 2004 Basel Meetings Paper}.

\bibitem[{CFO Forum and CRO Forum}, 2010]{CFO-CRO-2010}
{CFO Forum and CRO Forum} (2010).
\newblock {QIS 5 Technical Specification; Risk-free interest rates}.
\newblock {\em {Technical report, CFO Forum and CRO Forum.}}

\bibitem[Chibane et~al., 2009]{Chibane2009}
Chibane, M., Selvaraj, J., and Sheldon, G. (2009).
\newblock Building curves on a good basis.
\newblock working paper.

\bibitem[Chiu et~al., 2008]{Chiu_el2008}
Chiu, N.-C., Fang, S.-C., Lavery, J.~E., Lin, J.-Y., and Wang, Y. (2008).
\newblock Approximating term structure of interest rates using cubic {$L_1$}
  splines.
\newblock {\em European Journal of Operational Research}, 184(3):990--1004.

\bibitem[Cont, 2006]{Cont2006}
Cont, R. (2006).
\newblock Model uncertainty and its impact on the pricing of derivative
  instruments.
\newblock {\em Mathematical finance}, 16(3):519--547.

\bibitem[Cressie, 1990]{cressie1990origins}
Cressie, N. (1990).
\newblock The origins of kriging.
\newblock {\em Mathematical geology}, 22(3):239--252.

\bibitem[Cressie, 1993]{Cressie1993Wiley}
Cressie, N.~A.~C. (1993).
\newblock {\em {Statistics for Spatial Data, revised edition}}.
\newblock John Wiley \& Sons, New York.

\bibitem[Davis and Hobson, 2007]{Davis2007}
Davis, M.~H. and Hobson, D.~G. (2007).
\newblock The range of traded option prices.
\newblock {\em Mathematical Finance}, 17(1):1--14.

\bibitem[de~Andr{\'e}s~S{\'a}nchez and G{\'o}mez, 2004]{deAndresGomez2004}
de~Andr{\'e}s~S{\'a}nchez, J. and G{\'o}mez, A.~T. (2004).
\newblock Estimating a fuzzy term structure of interest rates using fuzzy
  regression techniques.
\newblock {\em European Journal of Operational Research}, 154(3):804--818.

\bibitem[Deb{\'o}n et~al., 2010]{Debon2010327}
Deb{\'o}n, A., Mart{\'\i}nez-Ruiz, F., and Montes, F. (2010).
\newblock {A geostatistical approach for dynamic life tables: The effect of
  mortality on remaining lifetime and annuities}.
\newblock {\em Insurance: Mathematics and Economics}, 47(3):327 -- 336.

\bibitem[Derman, 1996]{Derman1996}
Derman, E. (April 1996).
\newblock {Model Risk}.
\newblock {\em Quantitative Strategies Research Notes, Goldman Sachs}.

\bibitem[Eberlein and Jacod, 1997]{Eberlein1997}
Eberlein, E. and Jacod, J. (1997).
\newblock On the range of options prices.
\newblock {\em Finance and Stochastics}, 1(2):131--140.

\bibitem[El~Karoui et~al., 1998]{ElKaroui1998}
El~Karoui, N., Jeanblanc-Picqu{\`e}, M., and Shreve, S.~E. (1998).
\newblock {Robustness of the Black and Scholes formula}.
\newblock {\em Mathematical finance}, 8(2):93--126.

\bibitem[Fengler and Hin, 2015]{Fengler:2015}
Fengler, M.~R. and Hin, L.-Y. (2015).
\newblock A simple and general approach to fitting the discount curve under
  no-arbitrage constraints.
\newblock {\em Finance Research Letters}, 15:78--84.

\bibitem[Fernandes and Artes, 2015]{FernandesArtes2015}
Fernandes, G.~B. and Artes, R. (2015).
\newblock Spatial dependence in credit risk and its improvement in credit
  scoring.
\newblock forthcoming in European Journal of Operational Research.

\bibitem[Fries, 2013]{Fries2013}
Fries, C.~P. (2013).
\newblock {Curves and term structure models: Definition, calibration and
  application of rate curves and term structure models}.
\newblock DZ Bank AG; LMU Munich, Department of Mathematics.

\bibitem[Fujii et~al., 2010]{Fujii2010}
Fujii, M., Shimada, Y., and Takahashi, A. (2010).
\newblock Collateral posting and choice of collateral currency.
\newblock {\em Working Paper}.

\bibitem[Gan and Sheldon~Lin, 2015]{Gan2015138}
Gan, G. and Sheldon~Lin, X. (2015).
\newblock {Valuation of large variable annuity portfolios under nested
  simulation: A functional data approach}.
\newblock {\em Insurance: Mathematics and Economics}, 62(0):138 -- 150.

\bibitem[Gilli et~al., 2010]{gilli2010calibrating}
Gilli, M., Gro{\ss}e, S., and Schumann, E. (2010).
\newblock Calibrating the {Nelson-Siegel-Svensson} model.
\newblock {\em COMISEF working papers series}, (31).

\bibitem[Golchi et~al., 2015]{doi:10.1137/140976741}
Golchi, S., Bingham, D.~R., Chipman, H., and Campbell, D.~A. (2015).
\newblock {Monotone Emulation of Computer Experiments}.
\newblock {\em SIAM/ASA Journal on Uncertainty Quantification}, 3(1):370--392.

\bibitem[Green and Figlewski, 1999]{Green1999}
Green, T.~C. and Figlewski, S. (1999).
\newblock Market risk and model risk for a financial institution writing
  options.
\newblock {\em The Journal of Finance}, 54(4):1465--1499.

\bibitem[Guojun, 2013]{Gan2013795}
Guojun, G. (2013).
\newblock Application of data clustering and machine learning in variable
  annuity valuation.
\newblock {\em Insurance: Mathematics and Economics}, 53(3):795 -- 801.

\bibitem[Hagan and West, 2006]{Hagan2006}
Hagan, P.~S. and West, G. (2006).
\newblock Interpolation methods for curve construction.
\newblock {\em Applied Mathematical Finance}, 13(2):89--129.

\bibitem[Henaff, 2010]{Henaff2010}
Henaff, P. (2010).
\newblock A normalized measure of model risk.
\newblock working paper.

\bibitem[Hull and White, 2013]{Hull2013}
Hull, J. and White, A. (2013).
\newblock {LIBOR vs. OIS: The derivatives discounting dilemma}.
\newblock {\em Journal of Investment Management}, 11(3):14--27.

\bibitem[Iwashita, 2013]{Iwashita2013}
Iwashita, Y. (2013).
\newblock {Piecewise Polynomial Interpolations}.
\newblock Open Gamma Technical Report.

\bibitem[Jones et~al., 1998]{jones1998}
Jones, D.~R., Schonlau, M., and Welch, W. (1998).
\newblock {Efficient Global Optimization of Expensive Black-Box Functions}.
\newblock {\em Journal of Global Optimization}, 13(4):455--492.

\bibitem[Kanevski et~al., 2008]{Kanevski2008}
Kanevski, M., Maignan, M., Pozdnoukhov, A., and Timonin, V. (2008).
\newblock Interest rates mapping.
\newblock {\em Physica A: Statistical Mechanics and its Applications},
  387(15):3897--3903.

\bibitem[Kenyon and Stamm, 2012]{Kenyon2012}
Kenyon, C. and Stamm, R. (2012).
\newblock {\em {Discounting, Libor, CVA and Funding: Interest Rate and Credit
  Pricing}}.
\newblock Palgrave Macmillan.

\bibitem[Kerkhof and Melenberg, 2004]{kerkhof2004backtesting}
Kerkhof, J. and Melenberg, B. (2004).
\newblock Backtesting for risk-based regulatory capital.
\newblock {\em Journal of Banking \& Finance}, 28(8):1845--1865.

\bibitem[Kimeldorf and Wahba, 1970]{KW:1970}
Kimeldorf, G.~S. and Wahba, G. (1970).
\newblock A correspondence between bayesian estimation on stochastic processes
  and smooting by splines.
\newblock {\em The Annals of Mathematical Statistics}, 41(2):495--502.

\bibitem[Kleijnen and Van~Beers, 2012]{kleijnen2012monotonicity}
Kleijnen, J.~P. and Van~Beers, W.~C. (2012).
\newblock Monotonicity-preserving bootstrapped kriging metamodels for expensive
  simulations.
\newblock {\em Journal of the Operational Research Society}, 64(5):708--717.

\bibitem[Krige, 1951]{krige1951statistical}
Krige, D. (1951).
\newblock {\em {A statistical approach to some mine valuation and allied
  problems on the Witwatersrand: By DG Krige}}.
\newblock PhD thesis, University of the Witwatersrand.

\bibitem[Laurini and Moura, 2010]{Laurini:2010}
Laurini, M.~P. and Moura, M. (2010).
\newblock Constrained smoothing b-splines for the term structure of interest
  rates.
\newblock {\em Insurance: Mathematics and Economics}, 46(2):339--350.

\bibitem[Laurini and Ohashi, 2015]{LauriniOhashi2015}
Laurini, M.~P. and Ohashi, A. (2015).
\newblock A noisy principal component analysis for forward rate curves.
\newblock {\em European Journal of Operational Research}, 246(1):140--153.

\bibitem[Liu and Staum, 2010]{liu2010stochastic}
Liu, M. and Staum, J. (2010).
\newblock Stochastic kriging for efficient nested simulation of expected
  shortfall.
\newblock {\em Journal of Risk}, 12(3):3.

\bibitem[Ludkovski, 2015]{Ludkovski:2015}
Ludkovski, M. (2015).
\newblock Kriging metamodels for {B}ermudan option pricing.
\newblock working paper.

\bibitem[Maatouk and Bay, 2014a]{MaatoukMCQMC2014}
Maatouk, H. and Bay, X. (2014a).
\newblock {A New Rejection Sampling Method for Truncated Multivariate Gaussian
  Random Variables Restricted to Convex Sets}.
\newblock To appear in Monte Carlo and Quasi-Monte Carlo Methods 2014.
  Springer-Verlag, Berlin 2016.

\bibitem[Maatouk and Bay, 2014b]{maatouk:hal-01096751}
Maatouk, H. and Bay, X. (2014b).
\newblock {Gaussian Process Emulators for Computer Experiments with Inequality
  Constraints}.
\newblock in revision.

\bibitem[Maatouk et~al., 2015]{Maatouk201538}
Maatouk, H., Roustant, O., and Richet, Y. (2015).
\newblock {Cross-Validation Estimations of Hyper-Parameters of Gaussian
  Processes with Inequality Constraints}.
\newblock {\em Procedia Environmental Sciences}, 27:38 -- 44.
\newblock Spatial Statistics conference 2015.

\bibitem[Mardia et~al., 1996]{MKGL:1996}
Mardia, K.~V., Kent, J.~T., Goodall, C.~R., and Little, J.~A. (1996).
\newblock Kriging and splines with derivative information.
\newblock {\em Biometrika}, 83(1):207--221.

\bibitem[Matheron, 1963]{matheron1963principles}
Matheron, G. (1963).
\newblock Principles of geostatistics.
\newblock {\em Economic geology}, 58(8):1246--1266.

\bibitem[Morini, 2011]{Morini2011}
Morini, M. (2011).
\newblock {\em {Understanding and Managing Model Risk: A practical guide for
  quants, traders and validators}}.
\newblock Wiley.

\bibitem[Nelson and Siegel, 1987]{nelson1987parsimonious}
Nelson, C.~R. and Siegel, A.~F. (1987).
\newblock Parsimonious modeling of yield curves.
\newblock {\em The Journal of Business}, 60(4):pp. 473--489.

\bibitem[Pallavicini and Tarenghi, 2010]{pallavicini2010interest}
Pallavicini, A. and Tarenghi, M. (2010).
\newblock Interest-rate modeling with multiple yield curves.
\newblock {\em Available at SSRN 1629688}.

\bibitem[Paulson and Hart, 2006]{paulson2006spatial}
Paulson, N.~D. and Hart, C.~E. (2006).
\newblock {A spatial approach to addressing weather derivative basis risk: A
  drought insurance example}.
\newblock In {\em 2006 Annual Meeting of American Agricultural Economics
  Association}.

\bibitem[Ramponi, 2003]{ramponi2003adaptive}
Ramponi, A. (2003).
\newblock Adaptive and monotone spline estimation of the cross-sectional term
  structure.
\newblock {\em International Journal of Theoretical and Applied Finance},
  6(02):195--212.

\bibitem[Rasmussen and Williams, 2005]{Rasmussen:2005:GPM:1162254}
Rasmussen, C.~E. and Williams, C. K.~I. (2005).
\newblock {\em {Gaussian Processes for Machine Learning (Adaptive Computation
  and Machine Learning)}}.
\newblock The MIT Press.

\bibitem[Robert, 1995]{Robert}
Robert, C.~P. (1995).
\newblock Simulation of truncated normal variables.
\newblock {\em Statistics and Computing}, 5(2).

\bibitem[Roustant et~al.,
  2012]{Roustant:Ginsbourger:Deville:2012:JSSOBK:v51i01}
Roustant, O., Ginsbourger, D., and Deville, Y. (2012).
\newblock {DiceKriging, DiceOptim: Two R Packages for the Analysis of Computer
  Experiments by Kriging-Based Metamodeling and Optimization}.
\newblock {\em Journal of Statistical Software}, 51(1):1--55.

\bibitem[Santner et~al., 2003]{sant:will:notz:2003}
Santner, T.~J., Williams, B., and Notz, W. (2003).
\newblock {\em {The Design and Analysis of Computer Experiments}}.
\newblock Springer-Verlag.

\bibitem[Smith and Wilson, 2001]{SW:2001}
Smith, A. and Wilson, T. (2001).
\newblock Fitting yield curves with long term constraints.
\newblock Technical report, Technical report, Bacon and Woodrow.

\bibitem[Sousa et~al., 2012]{sousa2012machine}
Sousa, J.~B., Esqu{\'\i}vel, M.~L., and Gaspar, R.~M. (2012).
\newblock {Machine learning Vasicek model calibration with Gaussian processes}.
\newblock {\em Communications in Statistics-Simulation and Computation},
  41(6):776--786.

\bibitem[Steeley, 2008]{Steeley:2008}
Steeley, J.~M. (2008).
\newblock Testing term structure estimation methods: Evidence from the {UK
  STRIPS} market.
\newblock {\em Journal of Money, Credit and Banking}, 40(7):1489--1512.

\bibitem[Stutvoet, 2007]{stutvoet2007fitting}
Stutvoet, E. (2007).
\newblock {Fitting Financial Models to Market Data Using Kriging}.
\newblock Master's thesis, Delft university of Technology.

\bibitem[Svensson, 1994]{Svensson:1994}
Svensson, L.~E. (1994).
\newblock Estimating and interpreting forward interest rates: Sweden 1992-1994.
\newblock Technical report, National Bureau of Economic Research.

\end{thebibliography}
