A Note on Combining Machine Learning with Statistical Modeling for Financial Data Analysis

Abstract

This note revisits the ideas of the so-called semiparametric methods that we consider to be very useful when applying machine learning in insurance. To this aim, we first recall the main essence of semiparametrics like the mixing of global and local estimation and the combining of explicit modeling with purely data adaptive inference. Then, we discuss stepwise approaches with different ways of integrating machine learning. Furthermore, for the modeling of prior knowledge, we introduce classes of distribution families for financial data. The proposed procedures are illustrated with data on stock returns for five companies of the Spanish value-weighted index IBEX35.

1. Introduction

The Editors of this Special Issue pointed out that machine learning (ML) has no unanimous definition. In fact, the term “ML”, coined by Samuel (1959), is quite differently understood in the different communities. The general definition is that ML is concerned with the development of algorithms and techniques that allow computers to learn. The latter means a process of recognizing patterns in the data that are used to construct models; cf. “data mining” (Friedman 1998). These models are typically used for prediction. In this note, we speak about ML for data based prediction and/or estimation. In such a context, one may say that ML refers to algorithms that computer codes apply to perform estimation, prediction, or classification. As said, they rely on pattern recognition (Bishop 2006) for constructing purely data-driven models. The meaning of “model” is quite different here from what “model” or “modeling” means in the classic statistics literature; see Breiman (2001). He speaks of the data modeling culture (classic statistics) versus the algorithmic modeling coming from engineering and computer science. Many statisticians have been trying to reconcile these modeling paradigms; see Hastie et al. (2009). Even though the terminology comes from a different history and background, the outcome of this falls into the class of so-called semi-parametric methods; see Ruppert et al. (2003) or Härdle et al. (2004) for general reviews. According to that logic, ML would be a non-parametric estimation, whereas the explicit parametrization forms the modeling part from classic statistics.

Why should this be of interest? The Editors of this Special Issue also urge practitioners not to ignore what has already been learned about financial data when using presumably fully automatic ML methods. Regarding financial data for example, Buch-Larsen et al. (2005), Bolancé et al. (2012), Scholz et al. (2015, 2016), and Kyriakou et al. (2019) (among others) have shown the significant gains in estimation and prediction when including prior knowledge in nonparametric prediction. The first two showed how knowledge-driven data transformation improves nonparametric estimation of distribution and operational risk; the third paper used parametrically-guided ML for stock return prediction; the fourth imputed bond returns to improve stock return predictions; and the last proposed comparing different theory-driven benchmark models regarding their predictability. Grammig et al. (2020) combined the opposing modeling philosophies to predict stock risk premia. In this spirit, we will discuss the combination of purely data-driven methods with smart modeling, i.e., using prior knowledge. This will be exemplified along the analysis of the distributions (conditional and unconditional) of daily stock returns and calculations of their value-at-risk (VaR).

Notice finally that in the case of estimation, methods are desirable that permit practitioners to understand, maybe not perfectly, but quite well, what the method is doing to the data. This will facilitate the interpretation of results and any further inference. Admittedly, this is not always necessary and certainly also depends on the knowledge or imagination of the user. Yet, we believe it is often preferable to analyze data in a glass box than in a black box. This aspect is respected in our considerations.

Section 2 revisits the ideas of semiparametric statistics. Section 3 provides an intensive treatment of the distribution modeling followed by its combination with local smoothing. In Section 4, we give empirical illustrations. Section 5 concludes. In the Appendix are given additional details.

2. Preliminary Considerations and General Ideas

It is helpful to first distinguish between global and local estimation. Global means that the parameter or function applies at any point and to the whole sample. Local estimation applies only to a given neighborhood, like for kernel regression. It is clear that localizing renders a method much more flexible; however, the global part allows for an easy modeling, and its estimation can draw on the entire sample. Non-parametric estimators are not local by nature; for example, power series based estimators are not. Unless you want to estimate a function with discontinuities, local estimators are usually smoothing methods; see Härdle et al. (2004). This distinction holds also for complex methods (cf. the discussion about extensions of neural networks to those that recognize local features), which often turn out to be related to weighted nearest neighbor methods; see Lin and Jeon (2006) for random forests or Silverman (1984) for splines. The latter is already a situation where we face a mixture of global and local smoothing; another one is orthogonal wavelet series (Härdle et al. 1998).

Those mixtures are interesting because the global parts can borrow the strength from a larger sample and have a smoothing effect, while the local parts allow for the desired flexibility to detect local features. Power series can offer this only by including a (in practice unacceptable) huge number of parameters. This is actually a major problem of many complex methods, but mixtures allow substantially reducing this number. At the same time, they allow us to include prior knowledge about general features. For example, imposing shape restrictions is much simpler for mixtures (like splines) than it is for purely local smoothers; see Meyer (2008) and the references therein. Unless the number of parameters is pre-fixed, their selection happens via reduction through regularization, which can be implemented in many ways. Penalization methods like P-splines (Eilers et al. 2015) or LASSO (Tibshirani 1996) are popular. The corresponding problem for kernel, nearest neighbors, and related methods is the choice of the neighborhood size. In any case, one has to decide about the penalization criterion and a tuning parameter. The latter is until today an open question; presently, cross-validation-type methods are the most popular ones. For kernel based methods, see Heidenreich et al. (2013) and Köhler et al. (2014) for a review or Nielsen and Sperlich (2003) in the context of forecasting in finance.

The first question concerns the kind of prior information available, e.g., whether it is about the set of covariates, how they enter (linearly, additively, with interactions), the shape (skewness, monotonicity, number of modes, fat tails), or more generally, about smoothness. This is immediately followed by the question of how this can be included; in some cases, this is obvious (like if knowing the set of variables to be included); in some others, it is more involved (like including parameter information via Bayesian modeling). Knowledge about smoothness is typically supposed in order to justify a particular estimator and/or the selection method for the smoothing parameter. Information about the shape or how covariates enter the model comprises the typical ingredients of semiparametric modeling (Horowitz 1998) to improve nonparametric estimation (Glad 1998).

Consider the problem of estimating a distribution, starting with the unconditional case. In many situations, you are more interested in those regions for which data are hardly available. If you used then a standard local density estimator, you would try to estimate interesting parameters like VaR from only very few observations, maybe one to five, which is obviously not a good idea. Buch-Larsen et al. (2005) proposed to apply a parametric transformation using prior knowledge. Combining this way a local (kernel density) estimator with such a global one, however, allowed them to borrow strength from the model and from data that were further away. Similarly, consider conditional distributions. Locally, around a given value of the conditioning variable, you may have too few observations to estimate a distribution nonparametrically. Then, you may impose on this neighborhood the same probability law up to some moments, as we will do in our example below.1

Certainly, a good mixture of global and local fitting is problem-adapted. Then, the question falls into two parts: which is the appropriate parametric modeling, and how to integrate it with the flexible local estimator. For the former, you have to resort to expertise in the particular field. For the latter, we will discuss some popular approaches. All this will be exemplified with the challenges of modeling stock returns for five big Spanish companies and predicting their VaR.

In our example, the first step is to construct a parametric guide for the distribution of stock returns Y. To this aim, we introduce the class of generalized beta-generated (BG) distributions (going back, among others, to Eugene et al. (2002) and Jones (2004)), as this distribution class allows modeling skewed distributions with potentially long or fat tails. While this is not a completely new approach, we present it with an explicit focus on the above outlined objectives including the calculation of VaRs and combining it with nonparametric estimation and/or validation. Our validation is more related to model selection and testing, today well understood and established, and therefore kept short. The former, i.e., the combination with nonparametric estimation, is discussed for the problem of analyzing conditional distributions and can be extended to the combination with methods for estimating in high dimensions. For this example, in which the prior knowledge enters via a distribution class, we discuss two approaches: one is based on the method of moments, the other one on maximum likelihood. The latter is popular due to Rigby and Stasinopoulos (2005) and Severini and Staniswalis (1994). Rigby and Stasinopoulos (2005) considered a fully parametrized model in which each distribution parameter (potentially transformed with a known link) is written as an additive function of covariates, typically including a series of random effects. They proposed a backfitting algorithm (implemented in the R library GAMLSS) to maximize a penalized likelihood corresponding to a posterior mode estimation using empirical Bayesian arguments. Severini and Staniswalis (1994) started out with the parametric likelihood, but localized it by kernels. This is maximized then for some given values of the covariates.

3. A Practical Example

We now discuss the technical steps for the announced practical example. While this section focuses on the technical part, the empirical exercise is done in the next section.

3.1. Distribution Modeling

Often, the Student t distribution was used in financial econometrics and risk management to model the conditional asset returns, going back to Bollerslev (1987). However, it is well known that it does not describe very well the empirical features of most financial data. Therefore, several proposals have been made of skewed Student t distributions; see Theodossiou (1998) and Zhu and Galbraith (2010) for the context of finance or Jones and Faddy (2003) and Azzalini and Capitanio (2003) in (applied) statistics. For a more general discussion and compendium, see Rigby et al. (2019).

Let us consider two classes of skewed t distributions. Both are derived from the (generalized) BG classes; see Appendix A. These allow for the generation of many flexible distribution classes. The fist version depends on two parameters, whereas the second one depends on three. One may directly go for the second one. Remember, however, later on, that this is just a parametric guide for the semiparametric estimator. There is certainly a trade-off between the gain of flexibility, on the one hand, and the loss of its regularization, on the other hand. Moreover, each additional parameter in the global part may raise the costs of implementation and computation to an unacceptable degree. Both skewed t distributions are generated by taking a standard Student t as the baseline distribution in (A1) and (A2), respectively. Specifically, plugging in $F_Y(y;a,b)=\frac{1}{2}(1+y/\sqrt{a+b+y^2})$, a so-called scaled Student $t_2$, into (A1), gives our skewed t of Type 1 with the probability density function (pdf):

$$g_{T_1}(y;a,b)=k_1{(1+\frac{y}{\sqrt{a+b+y^2}})}^{a+1/2}{(1-\frac{y}{\sqrt{a+b+y^2}})}^{b+1/2},$$

(1)

where $k_1^{-1}=B(a,b)\sqrt{a+b}{2}^{a+b-1}$ and $y\in \mathbb{R}$. We write $Y\sim T_1(a,b)$; see Jones and Faddy (2003). Plugging the baseline cumulative density (cdf) into (A2) gives our skewed t of Type 2 with the pdf:

$$g_{T_2}(y;a,b,c)=k_2\frac{1}{{(a+b+y^2)}^{3/2}}{(1+\frac{y}{\sqrt{a+b+y^2}})}^{ac-1}{(1-\frac{1}{2^c}{(1+\frac{y}{\sqrt{a+b+y^2}})}^c)}^{b-1},$$

(2)

where $k_2=\frac{c(a+b)}{B(a,b)2^{ac}}$ and $y\in \mathbb{R}$. We write $Y\sim T_2(a,b,c)$; see Alexander and Sarabia (2010) and Alexander et al. (2012). The densities $g_{T_1}$, $g_{T_2}$ are illustrated in Figure 1.

For implementation, estimation, model selection, and in particular, for the possible integration of ML, we take a closer look at the basic properties of these two distributions. Note that for $a=b$ in (1), one obtains the classical Student t distribution with $2a$ degrees of freedom. The same is true for the three-parameter case (2) when setting $a=b$ with $c=1$. Furthermore, the cdfs are given by:

$$G_{T_1}(y;a,b)=I\left(\frac{1}{2}\right\{1+\frac{y}{\sqrt{a+b+y^2}}\};a,b),$$

(3)

and:

$$G_{T_2}(y;a,b,c)=I(F_Y^c(y;a,b);a,b),$$

(4)

respectively. Having explicit expressions for the pdf, one may think that likelihood based methods are straight-forward, including those for tests and model selection. Note, however, that for a sample of size m, you need to know the joint distribution of $(y_1,\dots ,y_m)$, which are not independent. In practice, many work with $log{\prod }_{i=1}^m{g}_{T_k}(y_i;\theta )$, ($k=1,2$ with $\theta =(a,b)$ for $k=1$, $\theta =(a,b,c)$ else) as a likelihood approximate, ignoring potential dependencies. This weakens both the efficiency of estimation and the validity of likelihood based inference. Therefore, it is interesting to look at alternatives. For estimation, recall the moment based approach; then, we need to to express parameters $a,b$, and if applicable c, in terms of estimable moments. For the case of the skewed t distribution of the first kind, we have:

$$E\left(Y^r\right)=\frac{{(a+b)}^{r/2}}{B(a,b)}\sum _{j=0}^r\left(\genfrac{}{}{0pt}{}{r}{j}\right)2^{-j}{(-1)}^{j}B(a-\frac{r}{2}-j,b-\frac{r}{2})$$

(5)

if $a,b>r/2$. An advantage of this expression is that we do not need to estimate the centered, standardized moments. This is even more an advantage when we try to estimate these moments nonparametrically. Similarly, for the skewed t distribution of the second kind:

$$E\left(Y^r\right)=\frac{{(a+b)}^{r/2}}{B(a,b)}\sum _{j=0}^r{(-1)}^j\left(\genfrac{}{}{0pt}{}{r}{j}\right)2^{-j}\sum _{k=0}^{\infty }\left(\genfrac{}{}{0pt}{}{-r/2}{k}\right){(-1)}^{k}B(a-\frac{r/2+j-k}{c},b).$$

(6)

Clearly, this reveals a problem because we cannot solve this (easily) for $a,b$, and c.

3.2. Financial Risk Measures

The parametric part provides us with explicit formulas that we can use to derive closed expressions for the risk measures, namely the VaR and tail moments. In fact, it can be shown that:

$$\begin{array}{ccc}\hfill {\mathrm{VaR}}_{T_1}[p;a,b]& =& \frac{\sqrt{a+b}(2{\mathrm{VaR}}_B[p;a,b]-1)}{2\sqrt{{\mathrm{VaR}}_B[p;a,b](1-{\mathrm{VaR}}_B[p;a,b])}},\mathrm{and}\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathrm{VaR}}_{T_2}[p;a,b,c]& =& \frac{\sqrt{a+b}(2{\mathrm{VaR}}_B^{1/c}[p;a,b]-1)}{2\sqrt{{\mathrm{VaR}}_B^{1/c}[p;a,b](1-{\mathrm{VaR}}_B^{1/c}[p;a,b])}},\hfill \end{array}$$

(8)

with $0\le p\le 1$, where ${\mathrm{VaR}}_B[p;a,b]$ denotes the VaR of a classical $\mathcal{B}e(a,b)$ distribution. Furthermore, if $Y\sim T_1(a,b)$, then the for a given $y_p$, its corresponding tail moments can we written as:

$$\begin{array}{ccc}\hfill E\left\{Y^k\right|Y\le y_p\}& =& E\{\frac{{(a+b)}^{k/2}{(2B-1)}^k}{2^k{B}^{k/2}{(1-B)}^{k/2}}|\frac{\sqrt{(a+b)}(2B-1)}{2\sqrt{B(1-B)}}\le y_p\}\hfill \\ & =& \frac{{(a+b)}^{k/2}}{2^k}\sum _{j=0}^k{a}_{j}E\{B^{k/2-j}{(1-B)}^{k/2}|B\le h\left(x_p\right)\},\hfill \end{array}$$

where $a_j={(-1)}^j\left(\genfrac{}{}{0pt}{}{k}{j}\right)2^{k-j}$ and $h\left(z\right)=\{a+b+z^2+\sqrt{(a+b)z^2+z^4}\}/\left\{2(a+b+z^2)\right\},$ which is increasing in z, with $a,b\in {\mathbb{R}}^{+}$. For $k=1$, we get the T(ail)VaR. If $B\sim \mathcal{B}e(a,b)$, then one has:

$$E\{B^{k/2-j}{(1-B)}^{k/2}|B\le h\left(x_p\right)\}=\frac{1-F_{\tilde{B}}\left(h\left(y_p\right)\right)}{F_B\left(h\left(y_p\right)\right)}\frac{B(k/2+a-j,b-k/2)}{B(a,b)},$$

(9)

where $\tilde{B}\sim \mathcal{B}e(k/2+a-j,b-k/2)$, $b>k/2$, and $h\left(z\right)$ as above. In sum, for the Type 1 class, we obtain explicit expressions for the tail moments that do not change when $(a,b)$ are functions of covariates.

3.3. Combining The Prior with Nonparametric Estimation

First, let us briefly consider the case of being only interested in the estimation of the unconditional distribution $g(·)$, say, in that of stock returns Y. Once a proper parametric choice (say $G_{\theta }$) is found and the parameters $\theta$ are estimated, we know that $\tilde{Y}:=G_{\widehat{\theta }}\left(Y\right)$ has a pretty smooth density $\tilde{g}$, which is not too far from the uniform $[0,1]$ distribution. Then, as $g\left(y\right)=\tilde{g}\left\{G_{\widehat{\theta }}\left(y\right)\right\}$, you obtain the final estimate. A next step to proceed could be to apply a nonparametric estimator like, for example, a kernel density for $\tilde{g}$. The prior estimation served here to stretch data where we had many observations, and contract them where we had only a few. You may say that this could also be done either by taking the empirical cdf for transformation or by taking the local bandwidth that always includes the same number of neighbors. This would be purely nonparametric and therefore renounce the use of prior information. In practice, these alternatives suffer from various problems like giving wiggly results, a much harder bandwidth choice, etc. For details and applications in actuarial sciences, consult Bolancé et al. (2012) and Martínez Miranda et al. (2009).

We turn now to the slightly more challenging problem of considering conditional distributions. For example, what is the stock return distribution of Telefónica and its VaR given a certain value of the Spanish IBEX35? Having huge datasets, you may try to estimate this with a fully nonparametric estimator. However, if you doubt that stationarity holds over a long period, you may want to restrict your dataset to not include more than twelve months (as an example). In such a case, you better resort to a parametric guide like above and turn $a,b$ (and c if applicable) into flexible functions of the conditioning covariate X. Again, a likelihood based approach may now look at $logpro{d}_{i=1}^m{g}_{T_k}(y_i;\theta \left(x_i\right))$ ($k=1,2$ with $\theta \left(x_i\right)=(a\left(x_i\right),b\left(x_i\right))$ for $k=1$, etc.) ignoring potential dependencies. Rigby and Stasinopoulos (2005) specified the elements of $\theta$ as additive, parametrized functions of $x_i$ with some random coefficients, and maximized a penalized version of this. In contrast, Severini and Staniswalis (1994) did not parametrize the elements of $\theta$, but maximized (along $\theta$) the smoothed version, i.e., $su{m}_{i=1}^{m}log{g}_{T_k}(y_i;\theta \left(x\right))K_h(x-x_i)$ for any x. Here, $K_h\left(v\right)=h^{-1}K(v/h)$ for a kernel function $K(·)$ with bandwidth h. That is, they obtained for given x an estimate of the value that $\theta$ takes at x. Notice that the latter method is usually implemented for $\theta$ containing only one element (typically the mean).

As said before, an interesting alternative is to estimate nonparametrically the moments of Y and then derive the elements of $\theta$. Thanks to Formula (5), we can do this for the Type 1 skewed t distributions. More specifically, the algorithm would look as follows: For sample $(x_i,y_i)$, $i=1,2,\dots ,m$ with $y_i$ being the return of the stock and $x_i$ the conditioning variable (in our application, the market return) conduct:

Step 1: Estimate the conditional first two moments ${\mu }_1$ and ${\mu }_2$ by:

$$\begin{array}{c}{\mu }_1\left(x_i\right)=E\left(y\right|x_i), \mathrm{and} log{\mu }_2\left(x_i\right)=E\left(y^2\right|x_i),\hfill \end{array}$$

where the nonparametric functions can be estimated, e.g., by kernel regression, splines, etc.

Step 2: You now may either take a grid over the range of X, with M grid points $x_j$, or you may calculate the estimates for each observation $x_i$. Let us call them ${\widehat{\mu }}_{1j}$ and ${\widehat{\mu }}_{2j}$, $j=1,2,\dots ,M$ for either case.

Step 3: Calculate estimates $({\widehat{a}}_j,{\widehat{b}}_j)$ by solving in $(a_j,b_j)$, $j=1,2,\dots ,M$, the non-linear system:

$$\begin{array}{ccc}\hfill {\widehat{\mu }}_{1j}& =& \frac{(a_j-b_j)\sqrt{a_j+b_j}}{2}·\frac{\Gamma (a_j-\frac{1}{2})\Gamma (b_j-\frac{1}{2})}{\Gamma \left(a_j\right)\Gamma \left(b_j\right)},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill {\widehat{\mu }}_{2j}& =& \frac{a_j+b_j}{4}·\frac{{(a_j-b_j)}^2+a_j-1+b_j-1}{(a_j-1)(b_j-1)}.\hfill \end{array}$$

(11)

Step 4: Then, the estimate of the conditional distribution is of the form:

$$\widehat{g}\left(y\right|x_j)=g_F\left(y\right|{\widehat{a}}_j,{\widehat{b}}_j).$$

Following Nielsen and Sperlich (2003), Kyriakou et al. (2019), and Mammen et al. (2019), you could use local linear regression in Step 1 for both functions, combined with the validated $R^2$ for the bandwidth choice. Obviously, any method known as ML, including LASSO variable selection, can be applied in this step. However, it is less evident how these methods could be combined with the aforementioned likelihood based approaches. Note further that the conditional distribution obtained by this strategy can also be used for semiparametric prediction of the unconditional distribution:

$$\widehat{g}\left(y\right)=\frac{1}{m}\sum _{i=1}^m{g}_F\left(y\right|{\widehat{a}}_i,{\widehat{b}}_i).$$

(12)

This shows that with ${\widehat{\mu }}_1$, ${\widehat{\mu }}_2$, you can predict the marginal distribution of Y for scenarios in which the distribution of X is changing (Dai et al. 2016). For example, we could predict the unconditional distribution of stocks for different distributions of the IBEX35. Note finally that Step 3 cannot easily be applied to the skewed t distribution of Type 2; recall (6). In such a case, you could only try to apply the idea of Rigby and Stasinopoulos (2005), but with procedures and algorithms that are still to be developed.

4. Empirical Illustration

Consider daily stock returns from 1 January 2015 to 31.12.2015 for five companies of the Spanish value-weighted index IBEX35; namely Amadeus (IT solutions for tourist industry), BBVA (global financial services), Mapfre (insurance market), Repsol (energy sector), and Telefónica (information and communications technology services); see

Table 1. Summary statistics. The sample size is $m=261$ for all sets.

Stocks	Amadeus	BBVA	Mapfre	Repsol	Telefónica
Maximum daily return	0.046286	0.040975	0.050847	0.073466	0.062264
Minimum daily return	−0.097367	−0.060703	−0.067901	−0.0877323	−0.051563
Mean	0.000900	−0.000452	−0.000623	−0.001416	−0.000408
Standard deviation	0.014601	0.016249	0.015942	0.021349	0.016301
Skewness	−1.163797	−0.465779	−0.723655	−0.166165	0.130372
Kurtosis	10.292160	3.824688	4.873980	5.435928	4.422885

. The returns are negatively skewed in four of the five companies considered, but positively skewed in the last one.

We first fit the data by the maximum likelihood (ignoring dependence) to both distribution classes, working with standardized data.

Table 2. Maximum likelihood estimates for the skewed t model of Type 1, standardized data. Standard errors are in parenthesis.

Stocks	Amadeus	BBVA	Mapfre	Repsol	Telefónica
$\widehat{a}$	6.194309	10.773980	7.271484	5.009976	7.083988
	(2.378890)	(8.474473)	(3.684818)	(1.980271)	(3.810294)
$\widehat{b}$	6.171897	10.76088	7.250156	5.005015	7.086958
	(2.378415)	(8.477441)	(3.686769)	(1.980433)	(3.810390)

and

Table 3. Maximum likelihood estimates for the skewed t model of Type 2, standardized data. Standard errors are in parenthesis.

Stocks	Amadeus	BBVA	Mapfre	Repsol	Telefónica
$\widehat{a}$	1.050617	0.935678	0.8685684	0.808572	1.120804
	(0.443072)	(0.407211)	(0.329048)	(0.363157)	(0.650834)
$\widehat{b}$	5.126098	7.144007	6.217545	2.998354	3.497796
	(2.091549)	(5.104088)	(3.184219)	(0.917090)	(1.110353)
$\widehat{c}$	2.973896	3.653026	3.617017	2.721519	2.331579
	(0.761879)	(0.733523)	(0.668503)	(0.698227)	(0.774010)

show the parameter estimates with standard errors.

To choose between these two models, one may use the Bayesian information criterion (BIC). However, as our (working) likelihood neglects the dependence structure, it might not be reliable. While the three-parameter model presents the largest values for BIC (not shown), the gain, however, is always close to, or smaller than, 1%.

An alternative is to apply ML methods comparing the parametric estimates with purely nonparametric ones. This is not recommendable any longer when switching to conditional distributions due to the moderate sample sizes. In Figure 2, you see how our models (in red) adapted to the empirical cdf (blue) for the stocks of Amadeus and BBVA. As expected, the three-parameter model gave slightly better fits. In practice, the interesting thing to see is where improvements occurred, if any. The practitioner has to judge then what is of interest for his/her problem; ML cannot do this for him/her. However, ML can offer specification tests; see Gonzales-Manteiga and Crujeiras (2013) for a review. For example, a test that formalizes our graphical analysis is the Kolmogorov–Smirnov ($KS$) test:

$$KS=\sup |F_m\left(y_i\right)-F(y_i;\widehat{\theta })|,i=1,2,\dots ,m,$$

where $F_m\left(y_i\right)$ is the empirical cdf and $F(y;\widehat{\theta })$ is the cdf of the particular model class with $\widehat{\theta }$ from Table 2 and Table 3. To calculate the p-value, we can use the parametric bootstrap:

Step 1: For the observed sample, find the maximum likelihood estimator, $F(y;\widehat{\mathit{\theta }})$, ${\widehat{F}}_m\left(y\right)$, and $KS$.

Step 2: Generate J bootstrap samples $y_1^{\left(j\right)},\dots ,y_m^{\left(j\right)}\sim F(y;\widehat{\mathit{\theta }})$ under $H_0$ (the data follow model F); fit them; and compute $F(y;{\widehat{\mathit{\theta }}}^{\left(j\right)})$, ${\widehat{F}}_m^{\left(j\right)}\left(y\right)$, and $K{S}^{\left(i\right)}$ for each bootstrap sample $j=1,2,\dots ,J$.

Step 3: Calculate the p-value as the fraction of synthetic bootstrap samples with a $KS$ statistic greater than the empirical $KS$ statistic obtained from the original data.

To obtain an approximate accuracy of the p-value for $ϵ=$ 0.01, we generated $J=\frac{1}{4}{ϵ}^{-2}=2500$ bootstrap samples.

Table 4. Bootstrap p-values for both models and all five datasets.

Stocks	Amadeus	BBVA	Mapfre	Repsol	Telefónica
Skewed t $T_1$	0.593	0.676	0.499	0.761	0.829
Skewed t $T_2$	0.732	0.908	0.733	0.732	0.915

shows the results for both distribution classes and all datasets. It can be seen that with all p-values larger than 0.499, both models could not be rejected at any reasonable significance level in any of the considered datasets.

Finally, let us see how different the VaR are, when calculated on the base of one model compared to the other; recall Equations (7) and (8). They were calculated at the 95% confidence level for all datasets; see

Table 5. Values at risk ${\mathrm{VaR}}_{T1}[0.05;a,b)]$ and ${\mathrm{VaR}}_{T2}[0.05;a,b,c]$ for the five stocks considered.

Stocks	Amadeus	BBVA	Mapfre	Repsol	Telefónica
$Va{R}_{T1}$	−0.024941	−0.028328	−0.028521	−0.040059	−0.029110
$Va{R}_{T2}$	−0.023089	−.02817	−0.027794	−0.038330	−0.028029

. The $T_1$ model with two parameters provided slightly higher VaR values than the $T_2$ model. The difference again seemed to be somewhat marginal, except maybe for Amadeus.

Let us turn to the estimation of the conditional distribution. Here, the integration of the ML happens by incorporating the covariates nonparametrically. For the sake of presentation and brevity, we restricted ourselves here to the moment based approach; for more details on the likelihood based one, we refer (besides the above cited literature) to the recent compendium of Rigby et al. (2019) and the references therein. Limiting ourselves to the moment based method automatically limited us to the skewed t1 class; recall Equation (6).2 Furthermore, for the sake of illustration, we limited the exercise to the estimation of all distribution parameters as nonparametric functions of one given covariate X, namely the IBEX35. It was obvious that estimates for ${\mu }_1$, ${\mu }_2$ could equally well be the result of a complex multivariate regression or a variable selection procedure like LASSO. We estimated ${\mu }_1$(IBEX35) and ${\mu }_2$(IBEX35) using different methods provided by standard software; the presented results were obtained from penalized (cubic) spline regression with data-driven penalization. For details, consult Ruppert et al. (2003) and the SemiPar project. Figure 3 gives an example of how this performed for the BBVA stocks.

Table 6. First two moments of stock returns for given IBEX35 values (looking at its quantiles).

IBEX35	Amadeus		BBVA		Mapfre		Repsol		Telefónica
Quartile	${\widehat{\mathit{\mu }}}_{1\mathit{j}}$	${\widehat{\mathit{\mu }}}_{2\mathit{j}}$	${\widehat{\mathit{\mu }}}_{1\mathit{j}}$	${\widehat{\mathit{\mu }}}_{2\mathit{j}}$	${\widehat{\mathit{\mu }}}_{1\mathit{j}}$	${\widehat{\mathit{\mu }}}_{2\mathit{j}}$	${\widehat{\mathit{\mu }}}_{2\mathit{j}}$	${\widehat{\mathit{\mu }}}_{1\mathit{j}}$	${\widehat{\mathit{\mu }}}_{2\mathit{j}}$	${\widehat{\mathit{\mu }}}_{1\mathit{j}}$
$Q_1$	−0.323500	2.336970	−0.493668	1.489340	−0.481228	2.215213	−0.430198	1.654110	−0.509551	1.941169
$Q_2$	0.025464	2.385443	0.094427	1.240531	0.060678	1.487417	0.056462	1.474321	0.034497	1.253373
$Q_3$	0.280358	2.492818	0.503731	1.523203	0.439138	1.623399	0.405799	1.698132	0.433284	1.421627

gives for the IBEX35 quantiles $Q_1=-0.00768$, $Median=Q_2=0.00079$, $Q_3=0.00687$ the corresponding moment estimates of the different stock returns;

Table 7. Parameter $(a_j,b_j)$ of the conditional stock return distributions for given IBEX35 values.

IBEX35	Amadeus		BBVA		Mapfre		Repsol		Telefónica
Quartile	${\widehat{\mathit{a}}}_{\mathit{j}}$	${\widehat{\mathit{b}}}_{\mathit{j}}$	${\widehat{\mathit{a}}}_{\mathit{j}}$	${\widehat{\mathit{b}}}_{\mathit{j}}$	${\widehat{\mathit{a}}}_{\mathit{j}}$	${\widehat{\mathit{b}}}_{\mathit{j}}$	${\widehat{\mathit{a}}}_{\mathit{j}}$	${\widehat{\mathit{b}}}_{\mathit{j}}$	${\widehat{\mathit{a}}}_{\mathit{j}}$	${\widehat{\mathit{b}}}_{\mathit{j}}$
$Q_1$	1.742906	2.136783	5.397247	6.904428	1.990996	2.690554	3.143186	4.048275	2.490487	3.398563
$Q_2$	1.736629	1.709150	5.492494	5.226327	3.133551	3.019128	3.184576	3.076590	5.016808	4.924298
$Q_3$	1.953575	1.639172	6.478997	5.008818	4.327494	3.356619	3.640170	2.851049	6.809153	5.483973

the corresponding $(a_j,b_j)$ for Formulas (10) and (11).

Table 8. The conditional value at risk $Va{R}_{T1}\left(IBEX35\right)$ for given IBEX35 values.

IBEX35	Amadeus	BBVA	Mapfre	Repsol	Telefónica
$Q_1$	−0.037728	−0.038910	−0.044802	−0.053801	−0.043939
$Q_2$	−0.031199	−0.027981	−0.030278	−0.041117	−0.029327
$Q_3$	−0.026029	−0.020886	−0.022744	−0.032601	−0.021913

gives the corresponding conditional VaRs obtained from Formula (7).

In Figure 4, you can see the entire conditional distributions for the three IBEX35 quantiles. Finally, in Figure 5, we plotted the resulting unconditional distributions when you integrate over the observed IBEX35 values; recall Equation (12). They reflect quite nicely the asymmetries and some fat tails.

5. Discussion and Conclusions

In this note, we revisited the ideas of semiparametric modeling to propose for ML what one could call glass box modeling. It integrates ML in a mixture of a global and a local part in which the global one is as a parametric guide for the nonparametric estimate. We discussed different advantages of such a kind of smart modeling and the steps to be performed. In our illustration (analyzing financial data), we proposed as a parametric guide some (generalized) beta-generated distributions. In particular, we considered two classes of skewed t distributions. This allowed us to work with analytical expressions for the pdf, cdf, moments, and quantile functions, including the VaR, even on a local level. An empirical application with five datasets of stock returns was performed for illustration.