Portfolio Optimization on Multivariate Regime-Switching GARCH Model with Normal Tempered Stable Innovation

Abstract

This paper uses simulation-based portfolio optimization to mitigate the left tail risk of the portfolio. The contribution is twofold. (i) We propose the Markov regime-switching GARCH model with multivariate normal tempered stable innovation (MRS-MNTS-GARCH) to accommodate fat tails, volatility clustering and regime switch. The volatility of each asset independently follows the regime-switch GARCH model, while the correlation of joint innovation of the GARCH models follows the Hidden Markov Model. (ii) We use tail risk measures, namely conditional value-at-risk (CVaR) and conditional drawdown-at-risk (CDaR), in the portfolio optimization. The optimization is performed with the sample paths simulated by the MRS-MNTS-GARCH model. We conduct an empirical study on the performance of optimal portfolios. Out-of-sample tests show that the optimal portfolios with tail measures outperform the optimal portfolio with standard deviation measure and the equally weighted portfolio in various performance measures. The out-of-sample performance of the optimal portfolios is also more robust to suboptimality on the efficient frontier.

1. Introduction

Empirical studies have found in the return of various financial instruments skewness and leptokurtotic—asymmetry, and higher peak around the mean with fat tails. Normal distribution has long been recognized as insufficient to accommodate these stylized facts, relying on which could drastically underestimate the tail risk of a portfolio. The $\alpha$-stable distribution is a candidate to incorporate this fact into decision-making, but its lack of finite moments could cause difficulties in modeling. The class of tempered stable distributions serves as a natural substitution, since it has finite moments while retaining many attractive properties of the $\alpha$-stable distribution. The class of tempered stable distributions is derived by tempering the tails of the $\alpha$-stable distribution (See Koponen 1995; Boyarchenko and Levendorskiĭ 2000; and Carr et al. 2002), or from the time changed Brownian motion. Since Barndorff-Nielsen and Levendorskii (2001) and Shephard et al. (2001) presented the normal tempered stable (NTS) distribution in finance using the time changed Brownian motion, it has been successfully applied for modeling the stock returns with high accuracy in literature including Kim et al. (2012), Kurosaki and Kim (2018), Anand et al. (2016) and Kim (2022).

Deviation from normality exists not only in raw returns, but also in filtered innovations of time series models. Since the foundational work in Engle (1982) and Bollerslev (1986), GARCH model has been widely applied to model volatility. This motivates us to use NTS distribution to accommodate the asymmetry, interdependence, and fat tails of the innovations of GARCH models. Several studies have found NTS distribution to be more favorable than other candidates. Kim et al. (2011) found that normal and $\mathit{t}$-distribution are rejected as innovation of GARCH model in hypothesis testing. Shi and Feng (2016) found that tempered stable distribution with exponential tempering yields better calibration of Markov Regime-Switching (MRS) GARCH model than $\mathit{t}$ or generalized error distribution. Generalized hyperbolic distribution, while very flexible, has too many parameters for estimation. NTS distribution has only three standardized parameters and is flexible enough to serve the purpose.

A drawback of GARCH model is the inadequate ability in modeling volatility spikes, which could indicate that the market is switching within regimes. Various types of MRS-GARCH models have been proposed, many of which suggest that the regime-switching GARCH model achieves a better fit for empirical data (for example, in Hamilton 1996). Marcucci (2005) finds that MRS-GARCH models outperform standard GARCH models in volatility forecasting in a time horizon shorter than one week. Naturally, the most flexible model allows all parameters to switch among regimes. However, as is shown in Henneke et al. (2011), the sampling procedure in the MCMC method for estimating such a model is time-consuming and renders it improper for practical use. We construct a new model based on the regime-switching GARCH model specified in Haas et al. (2004), which circumvents the path dependence problem in the Markov Chain model by specifying parallel GARCH models.

Bollerslev (1990) proposes the GARCH model with constant conditional correlation. The parsimonious DCC-GARCH model in Engle (2002) opens the door to modeling multivariate GARCH process with different persistence between variance and correlation. While flexible models allowing time-varying correlation are preferred in many cases, specification of multivariate GARCH models with both regime-specific correlations and variance dynamics involves a balance between flexibility and tractability. The univariate MRS-GARCH model in Haas et al. (2004) is generalized in Haas and Liu (2004) to a multivariate case. Unfortunately, it suffers from the curse of dimensionality, and thus is unsuitable for a high dimensional empirical study. In our model, we decompose variance and correlation so that the variance of each asset evolves independently according to a univariate MRS-GARCH model. The correlation is incorporated in the innovations modeled by a flexible Hidden Markov Model (HMM) that has multivariate NTS distribution as the conditional distribution in each regime. We will use the calibrated model to simulate returns for portfolio optimizations.

Modern portfolio theory is formulated as a trade-off between return and risk. The classical Markowitz Model finds the portfolio with the highest Sharpe ratio. However, variance has been criticized for not containing enough information on the tail of the distribution. Moreover, correlation is not sufficient to describe the interdependence of asset returns with non-elliptical distribution. Current regulations for financial business utilize the concept of Value at Risk (VaR), which is the percentile of the loss distribution, to model the risk of left tail events. There are several undesired properties that rendered it an insufficient criterion. First, it is not a coherent risk measure due to a lack of sub-additivity. Second, VaR is a non-convex and non-smooth function of positions with multiple local extrema for non-normal distributions, which causes difficulty in developing efficient optimization techniques. Third, a single percentile is insufficient to describe the tail behavior of a distribution, and thus might lead to an underestimation of risk.

Theory and algorithm for portfolio optimization with a conditional value-at-risk (CVaR) measure proposed in Krokhmal et al. (2001) and Rockafellar and Uryasev (2000) effectively addresses these problems. For continuous distributions, CVaR is defined by a conditional expectation of losses exceeding a certain VaR level. For discrete distributions, CVaR is defined by a weighted average of some VaR levels that exceed a specified VaR level. In this way, CVaR concerns both VaR and the losses exceeding VaR. As a convex function of asset allocation, a coherent risk measure and a more informative statistic, CVaR serves as an ideal alternative to VaR. A study on the comparison of the two measures can be found in Sarykalin et al. (2014).

Another closely related tail risk measure is conditional drawdown-at-risk (CDaR). Drawdown has been a constant concern for investors and is often used in the performance evaluation of a portfolio. It is much more difficult to climb out of a drawdown than drop into one, considering that it takes 100% return to recover from 50% relative drawdown. In behavioral finance, it is well documented that investors fear losses more than they value gains. However, the commonly used maximum drawdown only considers the worst scenario that only occur under some very special circumstances, and thus is very sensitive to the testing period and asset allocation. On the other hand, small drawdowns included in the calculation of average drawdown are acceptable and might be caused by pure noise, of which the minimization might result in overfitting. For instance, a Brownian motion would have drawdowns in a simulated sample path.

CDaR is proposed in Checkhlov et al. (2004) to address these concerns. While CVaR is determined by the distribution of return, CDaR is path-dependent. CDaR is essentially CVaR of the distribution of drawdowns. By this, we overcome the drawbacks of average drawdown and maximum drawdown. CDaR takes not only the depth of drawdowns into consideration, but also the length of them. Since the CDaR risk measure is the application of CVaR in a dynamic case, it naturally holds nice properties of CVaR, such as convexity with respect to asset allocation. The optimization method with constraints on CDaR has also been developed and studied in Checkhlov et al. (2004, 2005).

An optimization procedure could lead to an underperformed portfolio by magnifying the estimation error. Using historical sample paths as input implies an assumption that what happened in the past will happen in the future, which requires careful examination. It is found in Lim et al. (2011) that estimation errors of CVaR and the mean is large, resulting in unreliable allocation. An alternative way is to use multiple simulated returns. It is reasonable to expect outperformance during a crisis if we use the simulation of a model that captures extreme tail events as input for optimization with tail risk measures. Relevant research that resorts to copula does not address a regime switch (See Sahamkhadam et al. 2018).

We have identified several issues in GARCH-simulation-based portfolio optimization: deviation from normality, regime switch, high dimensional calibration and tail risk measures. This paper intends to address these issues in one shot. First, we propose the MRS-MNTS-GARCH model to accommodate fat tails, volatility clustering and regime switch. In the model, the volatility of each asset independently follows a regime-switch GARCH model, while the correlation of joint innovation of GARCH models follows a Hidden Markov Model. We specify the method for efficient model calibration. The MRS-MNTS-GARCH model is used to simulate sample paths. Then, the sample paths are used as input to portfolio optimization with tail risk measures CVaR and CDaR. We conduct in-sample tests sow show the goodness-of-fit. We conduct out-of-sample tests to show the outperformance of our approach in various measures. The optimal portfolios are also more robust to suboptimality of the optimization. The results suggest that, in practice, performance ratios can be improved compared to equally weighted portfolio by adopting simulation-based portfolio optimization with the MRS-MNTS-GARCH model. The empirical study suggests that using tail risk measures outperforms using standard deviation in terms of various performance measures and robustness.

The remainder of the paper is organized as follows. Section 2 introduces the preliminaries on NTS distribution and GARCH model. Section 3 specifies our model and methods for estimation and simulation. Section 4 introduces the portfolio optimization with tail risk measures. Section 5 is an empirical study on in-sample goodness-of-fit and out-of-sample performance in recent years.

2. Preliminaries

This section reviews the background knowledge on different risk measures, scenario-based estimation, normal tempered stable (NTS) distribution, and GARCH and regime-switching GARCH model.

2.1. Tail Risk Measures

In this subsection, we discuss and summarize popular tail risk measures in finance such as Value at Risk (VaR), Conditional VaR (CVaR), and conditional drawdown-at-risk (CDaR). Definitions of CDaR, and related properties of drawdown are following Checkhlov et al. (2004, 2005).

Consider a probability space $(\mathsf{\Omega },{\mathcal{F}}_{\infty },)$ and a portfolio with N assets. Suppose that $P_n$ is a stochastic process of a price for the n-th asset in the portfolio on the space

$$P_n:[0,\infty )\times \mathsf{\Omega }⟶\mathbb{R},n=1,2,\cdots ,N$$

with the real initial value $P_n(0,·)=P_{n,0}>0$.

The stochastic process of the compounded return of the n-th asset between time 0 to time $t>0$ is defined by

$${\mathcal{R}}_n(t,\omega )=\frac{P_n(t,\omega )-P_{n,0}}{P_{n,0}},t>0,$$

with ${\mathcal{R}}_n(0,·)=0$. Then, the compounded portfolio return $R(x,t,\omega )$ at time $t>0$ is equal to

$$\mathcal{R}(x,t,\omega )=\sum _{n=1}^N{x}_n{\mathcal{R}}_n(t,\omega ).$$

The definition of VaR for $\mathcal{R}(x,t,·)$ with the confidence level $\eta$ is

$${\mathrm{VaR}}_{\eta }\left(\mathcal{R}(x,t,·)\right)=min\left\{z\right|(\mathcal{R}(x,t,·)\le -z)\ge \eta \}.$$

If $\mathcal{R}(x,t,·)$ has a continuous cumulative distribution, then we have

$${\mathrm{VaR}}_{\eta }\left(\mathcal{R}(x,t,·)\right)=F_{-\mathcal{R}(x,t,·)}^{-1}(1-\eta ),$$

where $F_{-\mathcal{R}(x,t,·)}^{-1}$ is the inverse cumulative distribution function of $-\mathcal{R}(x,t,·)$. The definition of CVaR for $\mathcal{R}(x,t,·)$ with the confidence level $\eta$ is

$$\begin{array}{c}\hfill {\mathrm{CVaR}}_{\eta }\left(\mathcal{R}(x,t,·)\right)=\frac{1}{1-\eta }{\int }_{\eta }^1{\mathrm{VaR}}_z\left(\mathcal{R}(x,t,·)\right)dz.\end{array}$$

Equivalently, CVaR can be obtained by a minimization formula

$$\begin{array}{c}\hfill {\mathrm{CVaR}}_{\eta }\left(\mathcal{R}(x,t,·)\right)=\underset{\theta }{min}\left\{\theta +\frac{1}{\eta }E\left[max\{0,-\mathcal{R}(x,t,·)-\theta \}\right]\right\}.\end{array}$$

(1)

Following Checkhlov et al. (2004), denote the uncompounded return of the n-th asset by ${\mathcal{Q}}_n(t,w)$. The uncompounded portfolio return is defined by

$$\mathcal{Q}(x,t,\omega )=\sum _{n=1}^N{x}_n{\mathcal{Q}}_n(t,\omega ).$$

Let $T>0$, and $\omega \in \mathsf{\Omega }$, then ${\left(R(x,t,\omega )\right)}_{t\in [0,T]}=\left\{R(x,t,\omega )\right|t\in [0,T]\}$ is one sample path of portfolio return with the capital allocation vector x from time 0 to T. The drawdown ${\left(\mathrm{DD}(·,\omega )\right)}_{t\in [0,T]}$ of the portfolio return is defined by

$$\mathrm{DD}(t,\omega )=\underset{s\in [0,t]}{max}\mathcal{Q}(x,s,\omega )-\mathcal{Q}(x,t,\omega ),\mathrm{for}t\in [0,T]$$

(2)

is a risk measure assessing the decline from a historical peak in some variable, and $\left\{\mathrm{DD}\right(t,\omega \left)\right|t\in [0,T],\omega \in \mathsf{\Omega }\}$ is a stochastic process of the drawdown. The average drawdown (ADD) up to time T for $\omega \in \mathsf{\Omega }$ is defined by

$$\mathrm{ADD}(T,\omega )=\frac{1}{T}{\int }_0^T\mathrm{DD}(t,\omega )dt.$$

(3)

It is the time average of drawdowns observed between time 0 to T. The maximum drawdown (ADD) up to time T for $\omega \in \mathsf{\Omega }$ is defined by

$$\mathrm{MDD}(T,\omega )=\underset{t\in (0,T)}{max}\mathrm{DD}(t,\omega ).$$

(4)

It is the maximum of drawdowns that have occurred up to time T for $\omega \in \mathsf{\Omega }$.

Let

$$F_{\mathrm{DD}}(z,\omega )=\frac{1}{T}{\int }_0^T{1}_{\mathrm{DD}(t,\omega )\le z}dt$$

(5)

and

$${\zeta }_{\eta }\left(\omega \right)=\left\{\begin{array}{cc}inf\left\{z\right|F_{\mathrm{DD}}(z,\omega )\ge \eta \}\hfill & \mathrm{if}\eta \in (0,1]\hfill \\ 0\hfill & \mathrm{if}\eta =0\hfill \end{array}\right..$$

(6)

According to Checkhlov et al. (2005), the CDaR is defined as the CVaR of drawdowns. Applying (1) to drawdowns, we have

$${\mathrm{CDaR}}_{\eta }(T,\omega )=\underset{\zeta }{min}\left\{\zeta +\frac{1}{\eta T}{\int }_0^{T}max\left\{\mathrm{DD}(t,\omega )-\zeta ,0\right\}dt\right\}.$$

It can also be represented as

$$\begin{array}{c}{\mathrm{CDaR}}_{\eta }(T,\omega )\hfill \\ =\left(\frac{F_{\mathrm{DD}}({\zeta }_{\eta }\left(\omega \right),\omega )-\eta }{1-\eta }\right){\zeta }_{\eta }\left(\omega \right)+\frac{1}{(1-\eta )T}{\int }_0^T\mathrm{DD}(t,\omega )1_{\mathrm{DD}(t,\omega )>{\zeta }_{\eta }\left(\omega \right)}dt.\hfill \end{array}$$

(7)

Note that ${\mathrm{CDaR}}_0(T,\omega )={\mathrm{ADD}}_{\omega }\left(T\right)$ and ${\mathrm{CDaR}}_1(T,\omega )={\mathrm{MDD}}_{\omega }\left(T\right)$.

Let

$$F_{\mathrm{DD}}\left(z\right)=\frac{1}{T}E\left[{\int }_0^T{1}_{\mathrm{DD}(t,·)\le z}dt\right]$$

(8)

and

$${\zeta }_{\eta }=F_{\mathrm{DD}}^{-1}\left(\eta \right)=\left\{\begin{array}{cc}inf\left\{z\right|F_{\mathrm{DD}}\left(z\right)\ge \eta \}\hfill & \mathrm{if}\eta \in (0,1]\hfill \\ 0\hfill & \mathrm{if}\eta =0\hfill \end{array}\right..$$

(9)

Considering all the sample path $\left\{\mathcal{R}\right(x,t,\omega \left)\right|t\in [0,T],\omega \in \mathsf{\Omega }\}$, we define CDaR at $\eta$ by

$$\begin{array}{c}{\mathrm{CDaR}}_{\eta }\left(T\right)\hfill \\ =\left(\frac{F_{\mathrm{DD}}\left(\zeta \left(\eta \right)\right)-\eta }{1-\eta }\right)\zeta \left(\eta \right)+\frac{1}{(1-\eta )T}E\left[{\int }_0^T\mathrm{DD}(t,·)1_{\mathrm{DD}(t,·)>\zeta \left(\eta \right)}dt\right].\hfill \end{array}$$

(10)

2.2. Scenario-Based Estimation

In this subsection, we discuss the risk estimation using given scenarios (historical scenario or simulation). Suppose that the time interval $[0,T]$ is divided by a partition $\{0=t_0<t_1<\cdots <t_M=T\}$, and denote

$$R_m\left(x\right)=\mathcal{R}(x,t_m,·).$$

We select ${\omega }_s\in \mathsf{\Omega }$ where $s\in \{1,2,\cdots ,S\}$ and S is the number of scenarios. Then, we obtain scenarios of the portfolio return process $\mathcal{R}$ at time $t_m$:

$$R_m^s\left(x\right)=\mathcal{R}(x,t_m,{\omega }_s)$$

where $m\in \{0,1,2,\cdots ,M\}$, and $s\in \{1,2,\cdots ,S\}$. We calculate VaR, CVaR, DD, and CDaR under the simulated scenarios. The VaR with significant level $\eta$ under the simulated scenario is estimated by

$$\begin{array}{c}\hfill {\mathrm{VaR}}_{\eta }\left(R_m\left(x\right)\right)=-inf\left\{u|\frac{1}{S}\sum _{s=1}^S{1}_{R_m^s\left(x\right)<u}>1-\eta \right\}.\end{array}$$

(11)

Let $R_m^{\left(k\right)}\left(x\right)$ be the k-th smallest value in $\left\{R_m^s\left(x\right)\right|s=1,2,\cdots ,S\}$, then CVaR at the significant level $\eta$ is estimated according to the formula

$$\begin{array}{c}{\mathrm{CVaR}}_{\eta }\left(R_m\left(x\right)\right)\hfill \\ =-\frac{1}{\eta }\left(\frac{1}{S}\sum _{k=1}^{\lceil S\left(\eta \right)\rceil -1}{R}_m^{\left(k\right)}\left(x\right)+\left(\eta -\frac{\lceil S\eta \rceil -1}{S}\right)R_m^{(\lceil S\eta \rceil )}\left(x\right)\right),\hfill \end{array}$$

(12)

where $\lceil x\rceil$ denotes the smallest integer larger than x (See Rachev et al. 2008 for details).

The rate of return of the n-th asset between time $t_{i-1}$ and time $t_i$ is defined by

$$r_n(t_i,\omega )=\frac{P_n(t_i,\omega )-P_n(t_{i-1},\omega )}{P_n(t_{i-1},\omega )},t>0,$$

with $r_n(0,·)=0$. It is referred to as return in this paper. Let $x={(x_1,x_2,\cdots ,x_N)}^T\in {\mathbb{R}}^N$ be a capital allocation vector of a long-only portfolio satisfying ${\sum }_{n=1}^N{x}_n=1$ and $x_n\in [0,1]$ for all $n\in \{1,2,\cdots ,N\}$. Then, the portfolio return $U(x,t_i,{\omega }_s)$ between time $t_{i-1}$ and time $t_i$ is

$$U(x,t_i,{\omega }_s)=\sum _{n=1}^N{x}_n{r}_n(t_i,{\omega }_s).$$

The scenario-based estimation on uncompounded cumulative return at time $t_i$ is

$$Q(x,t_i,{\omega }_s)=\sum _{t=1}^t{U}_n(x,t_i,{\omega }_s).$$

Denote $Q(x,t_i,{\omega }_s)$ by $Q_i^s\left(x\right)$. Using (2)–(4), DD, ADD, and MDD on a single simulated scenario are estimated, respectively, by

$$\begin{array}{cc}\hfill {\mathrm{DD}}_{m,s}& :=\mathrm{DD}(t_m,{\omega }_s)=\underset{j\in \{0,1,2,\cdots ,m\}}{max}{Q}_j^s\left(x\right)-Q_m^s\left(x\right),\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mathrm{ADD}}_{M,s}& :=\mathrm{ADD}(t_M,{\omega }_s)=\frac{1}{M}\sum _{m=1}^M{\mathrm{DD}}_{m,s},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\mathrm{MDD}}_{M,s}& :=\mathrm{MDD}(t_M,{\omega }_s)=\underset{m\in \{1,2,\cdots ,M\}}{max}{\mathrm{DD}}_{m,s}.\hfill \end{array}$$

(15)

Furthermore, using (5)–(7), CDaR is estimated on a single simulated scenario by

$${\mathrm{CDaR}}_{\eta }(T,s):={\mathrm{CDaR}}_{\eta }(T,{\omega }_s)=\frac{1}{(1-\eta )M}\sum _{m=1}^M{\mathrm{DD}}_{m,s}{1}_{{\mathrm{DD}}_{m,s}>{\zeta }_{\eta }^s},$$

where

$${\zeta }_{\eta }^s:={\zeta }_{\eta }\left(w_s\right)=\left\{\begin{array}{cc}inf\left\{z\right|F_{\mathrm{DD}}(z,s)\ge \eta \}\hfill & \mathrm{if}\eta \in (0,1]\hfill \\ 0\hfill & \mathrm{if}\eta =0\hfill \end{array}\right..$$

with

$$F_{\mathrm{DD}}(z,s)=\frac{1}{M}\sum _{m=1}^M{1}_{{\mathrm{DD}}_{m,s}\le z}.$$

Equivalently, it is estimated by the optimization formula

$${\mathrm{CDaR}}_{\eta }(T,s)=\underset{\zeta }{min}\left\{\zeta +\frac{1}{(1-\eta )M}\sum _{m=1}^{M}max\left({\mathrm{DD}}_{m,s}-\zeta ,0\right)\right\}.$$

With (8)–(10), DD and CDaR on multiple scenarios are estimated by

$$\mathrm{DD}\left(m\right)=\frac{1}{S}\sum _{s=1}^S{\mathrm{DD}}_{m,s}$$

and

$$\begin{array}{c}{\mathrm{CDaR}}_{\eta }\left(T\right)\hfill \\ =\left(\frac{F_{\mathrm{DD}}\left(\zeta \left(\eta \right)\right)-\eta }{1-\eta }\right)\zeta \left(\eta \right)+\frac{1}{(1-\eta )M}\sum _{s=1}^S\sum _{m=1}^M\frac{1}{S}{\mathrm{DD}}_{m,s}{1}_{{\mathrm{DD}}_{m,s}>\zeta \left(\eta \right)},\hfill \end{array}$$

(16)

respectively, where

$$F_{\mathrm{DD}}\left(z\right)=\frac{1}{S}\sum _{s=1}^S{F}_{\mathrm{DD}}(z,s)$$

and

$${\zeta }_{\eta }=\left\{\begin{array}{cc}inf\left\{z|F_{\mathrm{DD}}\left(z\right)\ge \eta \right\}\hfill & \mathrm{if}\eta \in (0,1]\hfill \\ 0\hfill & \mathrm{if}\eta =0\hfill \end{array}\right..$$

In this paper, we will apply time series models to generate the scenario of $r_n(m,s)$ for $m\in \{1,2,\cdots ,M\}$, $s\in \{1,2,\cdots ,S\}$ and $n\in \{1,2,\cdots ,N\}$.

2.3. Multivariate Normal Tempered Stable Distribution

Let $\mathcal{T}$ be a strictly positive random variable defined by the characteristic function for $\lambda \in (0,2)$ and $\theta >0$

$${\varphi }_{\mathcal{T}}\left(u\right)=exp\left(-\frac{2{\theta }^{1-\frac{\lambda }{2}}}{\lambda }({(\theta -iu)}^{\frac{\lambda }{2}}-{\theta }^{\frac{\lambda }{2}})\right).$$

Then, $\mathcal{T}$ is referred to as a tempered stable subordinator with parameters $\lambda$ and $\theta$. Suppose that $\mathbf{\xi }={({\xi }_1,\cdots ,{\xi }_N)}^T$ is a N-dimensional standard normal distributed random vector with a $N\times N$ covariance matrix $\mathsf{\Sigma }$, i.e., $\mathbf{\xi }\sim \Phi (0,\mathsf{\Sigma })$. The N-dimensional NTS distributed random vector $X={(X_1,\cdots ,X_N)}^T$ is defined by

$$\mathbf{X}=\mathbf{\mu }+\mathbf{\nu }(T-1)+\sqrt{T}\mathrm{diag}\left(\mathbf{\gamma }\right)\mathbf{\xi },$$

where $\mathbf{\mu }={({\mu }_1,\cdots ,{\mu }_N)}^T\in {\mathbb{R}}^N$, $\mathbf{\nu }={({\nu }_1,\cdots ,{\nu }_N)}^T\in {\mathbb{R}}^N$, $\mathbf{\gamma }={({\gamma }_1,\cdots ,{\gamma }_N)}^T\in {\mathbb{R}}_{+}^N$ and $\mathcal{T}$ is independent of $\mathbf{\xi }$. The N-dimensional multivariate NTS distribution specified above is denoted by

$$X\sim {\mathrm{MNTS}}_N(\lambda ,\theta ,\mathbf{\nu },\mathbf{\gamma },\mathbf{\mu },\mathsf{\Sigma }).$$

Let ${\mu }_n=0$, and ${\gamma }_n=\sqrt{1-{\nu }_n^2\left(\frac{2-\lambda }{2\theta }\right)}$, $\left|{\nu }_n\right|<\sqrt{\frac{2\theta }{2-\lambda }}$ for $n\in \{1,\cdots ,N\}$. This yields a random variable having zero mean and unit variance. In this case, the random variable is referred to as standard NTS distributed random variable and denoted by $X\sim {\mathrm{stdMNTS}}_N(\lambda ,\theta ,\mathbf{\nu },\mathsf{\Sigma })$. The covariance matrix of X is given by

$${\mathsf{\Sigma }}_X=\mathrm{diag}\left(\mathbf{\gamma }\right)\mathsf{\Sigma }\mathrm{diag}\left(\mathbf{\gamma }\right)+\frac{2-\lambda }{2\theta }{\mathbf{\nu }}^{\top }\mathbf{\nu }.$$

2.4. Regime-Switching GARCH Model

GARCH($p,q$) model has been studied extensively as a model for volatility.

$$\begin{array}{cc}\hfill r_t& =\eta +{\sigma }_t{ϵ}_t,\hfill \\ \hfill u_t& ={\sigma }_t{ϵ}_t,\hfill \\ \hfill {\sigma }_t^2& =\omega +\sum _{i=1}^p{\alpha }_i{u}_{t-i}^2+\sum _{i=1}^q{\beta }_i{\sigma }_{t-i}^2.\hfill \end{array}$$

where $r_t$ is the return at time t, ${\sigma }_t$ is the variance at time t, $\eta ,\omega ,{\alpha }_i,{\beta }_i$ are parameters. We will use GARCH(1,1) in this paper.

Before specifying our model, it would be clear to first specify the univariate Markov regime-switching model in Haas et al. (2004) as follows

$$\begin{array}{cc}\hfill r_t& ={\eta }_{{\Delta }_t}+{\sigma }_{{\Delta }_t,t}{ϵ}_t,\hfill \\ \hfill u_t& ={\sigma }_{{\Delta }_t,t}{ϵ}_t,\hfill \\ \hfill {\mathbf{\sigma }}_t^2& =\mathbf{\omega }+\mathbf{\alpha }{u}_{t-1}^2+\mathbf{\beta }\circ {\mathbf{\sigma }}_{t-1}^2,\hfill \\ \hfill {ϵ}_t& \stackrel{\mathrm{iid}}{\sim }\mathrm{N}(0,1),\hfill \end{array}$$

(17)

where

• $k=$ number of regimes in a Markov chain;
• ${\Delta }_t=\{1,\cdots ,k\}=$ Markov chain with an irreducible and primitive transition matrix that determines the regime at time t;
• $r_t=$ return at time t;
• ${\mathbf{\sigma }}_t^2=({\sigma }_{1,t}^2,\cdots ,{\sigma }_{k,t}^2)=$ vector of parallel variance;
• $\mathbf{\omega }=({\omega }_1,\cdots ,{\omega }_k)$, $\mathbf{\alpha }=({\alpha }_1,\cdots ,{\alpha }_k)$, $\mathbf{\beta }=({\beta }_1,\cdots ,{\beta }_k)$. $\mathbf{\omega }$, $\mathbf{\alpha }$, $\mathbf{\beta }$ are vectors of coefficients;
• ${\eta }_{{\Delta }_t}=$ regime specific mean return in regime ${\Delta }_t$ at time t;
• ${\sigma }_{{\Delta }_t,t}=$ the standard deviation in regime ${\Delta }_t$ at time t;
• ∘ denotes element-wise product.

There are k parallel GARCH models evolving simultaneously. The Markov chain ${\Delta }_t$ determines which GARCH model is realized at next moment.

To better understand (17), let us create the following matrix

$${\left[{\sigma }_{{\Delta }_t,t}\right]}_{{\Delta }_t=1,\dots ,k,t=1,\dots T}=\left[\begin{array}{ccccc}{\sigma }_{1,1}& {\sigma }_{1,2}& \cdots & {\sigma }_{1,t}& \cdots \\ {\sigma }_{2,1}& {\sigma }_{2,2}& \cdots & {\sigma }_{1,t}& \cdots \\ ⋮& ⋮& \ddots & \cdots & \cdots \\ {\sigma }_{k,1}& {\sigma }_{k,2}& \cdots & {\sigma }_{k,t}& \cdots \end{array}\right].$$

Each column is a total of k parallel standard deviations at time t. Each row is the process of one regime through time 1 to T. The process of standard deviation is generated as follows:

Step 1: At $t=1$, the initial column is generated.

Step 2: ${\Delta }_1$ decides which row is realized.

For example, suppose ${\Delta }_1=2$. That is, ${\sigma }_{2,1}$ is realized in real-world at $t=1$.

Step 3: By vector calculation, 2nd column is generated.

For example, suppose ${\sigma }_{2,1}$ is realized at $t=1$. We calculate

$$\left[\begin{array}{c}{\sigma }_{1,2}\\ {\sigma }_{2,2}\\ ⋮\\ {\sigma }_{k,2}\end{array}\right]=\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_k\end{array}\right]+{\sigma }_{2,1}{ϵ}_1\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_k\end{array}\right]+\left[\begin{array}{c}{\beta }_1\\ {\beta }_2\\ ⋮\\ {\beta }_k\end{array}\right]\circ \left[\begin{array}{c}{\sigma }_{1,1}\\ {\sigma }_{2,1}\\ ⋮\\ {\sigma }_{k,1}\end{array}\right].$$

Step 4: ${\Delta }_2$ Decides which row is realized.

For example, suppose ${\Delta }_2=1$. That is, ${\sigma }_{1,2}$ is realized in real world at $t=2$. A regime switch from regime 2 to regime 1 takes place.

Step 5: Repeat Step 1 to Step 4.

It is important to note that in the example, the realized variance at time $t=2$, ${\sigma }_{1,2}$, is calculated with both the realized variance at time $t=1$, ${\sigma }_{2,1}$, and unrealized ${\sigma }_{1,1}$. This is obvious from the vector equation shown above. ${\sigma }_{1,1}$ is in the parallel process and is not realized in real world.

More generally, when there is a regime switch, from the vector equation ${\mathbf{\sigma }}_t^2=\mathbf{\omega }+\mathbf{\alpha }{u_t}^2+\mathbf{\beta }\circ {\mathbf{\sigma }}_{t-1}^2$, we know that the variance ${\sigma }_{{\Delta }_t,t}$ is determined by ${\sigma }_{{\Delta }_t,t-1}$ rather than ${\sigma }_{{\Delta }_{t-1},t-1}$. ${\sigma }_{{\Delta }_t,t-1}$ is a variance in regime ${\Delta }_t$ at time t, which is generated simultaneously with ${\sigma }_{{\Delta }_{t-1},t-1}$, but does not exist in reality at time t or $t-1$. This makes the model different from the one in Henneke et al. (2011) that requires all parameters to switch regime according to a Markov chain. In that model, there is no parallel process, and ${\sigma }_{{\Delta }_t,t}$ is determined by ${\sigma }_{{\Delta }_{t-1},t-1}$.

The stationary conditions derived in Haas et al. (2004) require definitions of the matrices:

$$M_{ji}=p_{ji}(\mathbf{\beta }+\alpha e_i^{\top }),i,j=1,\dots ,k$$

and the block matrix

$$M=\left[M_{ji}\right],i,j=1,\dots ,k.$$

The necessary and sufficient condition for stationarity is $\rho \left(M\right)<1$, where $\rho \left(M\right)$ denotes the largest eigenvalue of matrix M.

3. MRS-MNTS-GARCH Model

This section defines the MRS-MNTS-GARCH model and introduces the procedures for model fitting and sample path simulation.

3.1. Model Specification

To model a multivariate process, we assume that the variance of each individual asset evolves independently according to the univariate regime-switching GARCH model specified previously with possibly different number of regimes, while the correlation between assets is modeled separately by stdMNTS distribution in the joint standard innovations. The joint innovations is defined bya HMM with stdMNTS as conditional distribution in each regime.

The MRS-MNTS-GARCH model is defined by

$$\begin{array}{cc}\hfill {\mathbf{r}}_t& ={\mathbf{\eta }}_t+{\mathbf{\sigma }}_t\circ {\mathbf{ϵ}}_t\hfill \\ \hfill {\mathbf{u}}_t& ={\mathbf{\sigma }}_t\circ {\mathbf{ϵ}}_t\hfill \\ \hfill {\left({\mathbf{\sigma }}_t^{\left(i\right)}\right)}^2& ={\omega }^{\left(i\right)}+{\alpha }^{\left(i\right)}{\left(u_{t-1}^{\left(i\right)}\right)}^2+{\beta }^{\left(i\right)}\circ {\left({\mathbf{\sigma }}_{t-1}^{\left(i\right)}\right)}^2,i=1,\dots ,N\hfill \\ \hfill {\mathbf{ϵ}}_t& \stackrel{\mathrm{iid}}{\sim }\mathrm{stdMNTS}({\lambda }_{{\Delta }_t^D},{\theta }_{{\Delta }_t^D},{\mathbf{\nu }}_{{\Delta }_t^D},{\mathsf{\Sigma }}_{{\Delta }_t^D}),\hfill \end{array}$$

where

• $N=$ number of assets;
• $k^{\left(i\right)}=$ number of regimes of the i-th asset;
• ${\Delta }^{\left(i\right)}=\{1,\cdots ,k^{\left(i\right)}\}=$ Markov chain with an irreducible and primitive transition matrix that determines the regime of the i-th asset at time t;
• ${\Delta }_t^{\left(i\right)}=\{1,\cdots ,k^{\left(i\right)}\}=$ realization of ${\Delta }^{\left(i\right)}$ at time t;
• ${\Delta }^D=\{1,\cdots ,k^D\}=$ Markov chain with an irreducible and primitive transition matrix that determines the regime of the joint innovations at time t;
• ${\Delta }_t^D=\{1,\cdots ,k^D\}=$ realization of ${\Delta }^D$ at time t;
• ${\mathbf{r}}_t=(r_t^{\left(1\right)},\dots ,r_t^{\left(N\right)})=$ vector return of N assets at time t;
• ${\mathbf{\eta }}_t=({\eta }_{{\Delta }_t^{\left(1\right)}}^{\left(1\right)},\dots ,{\eta }_{{\Delta }_t^{\left(N\right)}}^{\left(N\right)})=$ vector mean return at time t;
• ${\eta }_{{\Delta }_t^{\left(i\right)}}^{\left(i\right)}=$ mean return of the i-th asset at time t, ${\eta }_{{\Delta }_t^{\left(i\right)}}^{\left(i\right)}$ is determined by the regime of the i-th asset at time t, ${\Delta }_t^{\left(i\right)}$;
• ${\mathbf{\sigma }}_t=\left({\sigma }_{{\Delta }_t^{\left(1\right)},t}^{\left(1\right)},\dots ,{\sigma }_{{\Delta }_t^{\left(N\right)},t}^{\left(N\right)}\right)=$ vector of standard deviation of N assets at time t;
• ${\mathbf{\sigma }}_t^{\left(i\right)}=\left({\sigma }_{1,t}^{\left(i\right)},\dots ,{\sigma }_{k^{\left(i\right)},t}^{\left(i\right)}\right)=$ standard deviation vector of the i-th asset with ${\left({\mathbf{\sigma }}_t^{\left(i\right)}\right)}^2=\left({\left({\sigma }_{1,t}^{\left(i\right)}\right)}^2,\dots ,{\left({\sigma }_{k^{\left(i\right)},t}^{\left(i\right)}\right)}^2\right)$;
• $u_{t-1}^{\left(i\right)}=$ the i-th element of ${\mathbf{u}}_t$, with ${\mathbf{u}}_t=\left(u_t^{\left(1\right)},\dots ,u_t^{\left(N\right)}\right)$;
• ${\mathbf{\omega }}^{\left(i\right)}=({\omega }_1^{\left(i\right)},\dots ,{\omega }_{k^{\left(i\right)}}^{\left(i\right)})$, ${\mathbf{\alpha }}^{\left(i\right)}=({\alpha }_1^{\left(i\right)},\dots ,{\alpha }_{k^{\left(i\right)}}^{\left(i\right)})$, ${\mathbf{\beta }}^{\left(i\right)}=({\beta }_1^{\left(i\right)},\dots ,{\beta }_{k^{\left(i\right)}}^{\left(i\right)})$. ${\mathbf{\omega }}^{\left(i\right)}$, ${\mathbf{\alpha }}^{\left(i\right)}$, ${\mathbf{\beta }}^{\left(i\right)}$ are the vector of coefficients of the i-th asset in $k^{\left(i\right)}$ regimes;
• ${\mathbf{ϵ}}_t=({ϵ}_t^{\left(1\right)},\dots ,{ϵ}_t^{\left(N\right)})=$ joint innovations;
• ∘ denotes element-wise product.

For each asset, the volatility process is the same as that of the univariate regime-switching GARCH model in (17). It is a classical GARCH process when it does not shift within regimes. Unlike the most flexible model in which the variance at time t is calculated with the variance at time $t-1$ in the last regime, when a regime shift takes place at time t, the variance in the model at time t is determined by the variance at time $t-1$ within the new regime. The regimes of each asset and joint innovations evolve independently, mitigating the difficulty in estimation caused by all parameters switching regimes simultaneously.

3.2. Parameter Estimation

Although NTS distribution is flexible enough to accommodate many stylized facts of asset return, the absence of an analytical density function presents some difficulties in estimation. Our estimation methodology is adapted from Kim et al. (2011, 2012), which are among the first works to incorporate NTS distribution in portfolio optimization. Another estimation of stdMNTS distribution with EM algorithm combined with fast Fourier transform is studied in Bianchi et al. (2016). Below we describe the procedure to fit the model.

First, we fit the univariate Markov regime-switching GARCH model with t innovation on the index return and extract the innovations and the regimes at each time t, ${\Delta }_t^D$. It is intuitive to assume that the correlation between assets change drastically when a market regime shift takes place. Thus, here we assume that the Markov chain ${\Delta }^D$ that determines the regime of innovations is the same as the Markov chain that determines the regime of the market index. That is, the sample path of the regime of the innovation is derived by extracting the sample path of the regime of the market index. The innovations of the model calibrated by the market index are used to estimate tail parameters ${\lambda }_{{\Delta }^D},{\theta }_{{\Delta }^D}$ in each regime.

Next, we fit the univariate model on each asset and extract the innovations for subsequent estimation of stdMNTS distribution. In each regime, the stdMNTS distribution has two tail parameters ${\lambda }_{{\Delta }^D},{\theta }_{{\Delta }^D}$ and one skewness parameter vector ${\mathbf{\nu }}_{{\Delta }^D}$. Common tail parameters ${\lambda }_{{\Delta }^D},{\theta }_{{\Delta }^D}$ are assumed for individual constituents, which are estimated from the market index in each regime. The skewness parameter ${\mathbf{\nu }}_{{\Delta }^D}$ is calibrated by curve fitting for each asset.

Finally, we plug the empirical covariance matrix of the innovations in each regime ${{\mathsf{\Sigma }}_X}_{{\Delta }^D}$ in the formula ${\mathsf{\Sigma }}_{{\Delta }^D}=\mathrm{diag}{\left({\mathbf{\gamma }}_{{\Delta }^D}\right)}^{-1}({{\mathsf{\Sigma }}_X}_{{\Delta }^D}-\frac{2-{\lambda }_{{\Delta }^D}}{2{\theta }_{{\Delta }^D}}{\mathbf{\nu }}_{{\Delta }^D}^{\top }{\mathbf{\nu }}_{{\Delta }^D})\mathrm{diag}{\left({\mathbf{\gamma }}_{{\Delta }^D}\right)}^{-1}$ to compute ${\mathsf{\Sigma }}_{{\Delta }^D}$. Note that ${\mathsf{\Sigma }}_{{\Delta }^D}$ is supposed to be a positive definite matrix. To guarantee this, we can either substitute ${\mathbf{\nu }}_{{\Delta }^D}$ with the closest vector to ${\mathbf{\nu }}_{{\Delta }^D}$ in $L^2$ norm which renders ${{\mathsf{\Sigma }}_X}_{{\Delta }^D}-\frac{2-{\lambda }_{{\Delta }^D}}{2{\theta }_{{\Delta }^D}}{\mathbf{\nu }}_{{\Delta }^D}^{\top }{\mathbf{\nu }}_{{\Delta }^D}$ positive definite matrix, or directly find the closest positive definite matrix to ${\mathsf{\Sigma }}_{{\Delta }^D}$. We choose the second method. Furthermore, the usual issue of estimating covariance matrix arises here as well, especially when a regime only has a short lasting period. To address this issue, we apply the denoising technique in Laloux et al. (1999) on the positive definite matrix derived earlier to estimate ${\mathsf{\Sigma }}_{{\Delta }^D}$.

The number of regimes of each asset is set as equal to that of the market index. We find that univariate Markov switching model with more than three regimes often has one unstable regime that lasts for a very short period and switches frequently and sometimes has one regime with clearly bimodal innovations distribution. Thus, it is desirable to limit the number of regimes smaller than 4, check unimodal distribution with Dip test and choose the one with the highest BIC value.

The procedure is summarized as follows.

Step 1: Fit the univariate Markov regime-switching GARCH model with t innovation on the index return and extract the innovations and the sample path of regimes.

Step 2: Estimate common tail parameters ${\lambda }_{{\Delta }^D},{\theta }_{{\Delta }^D}$ in each regime.

Step 3: Fit univariate model on each asset and extract the innovations.

Step 4: Estimate skewness parameters ${\nu }_{{\Delta }^D}$ for assets in each regime, calculate ${\gamma }_{{\Delta }^D}=\sqrt{1-{{\nu }_{{\Delta }^D}}^2\left(\frac{2-{\lambda }_{{\Delta }^D}}{2{\theta }_{{\Delta }^D}}\right)}$.

Step 5: Calculate the matrix $\mathrm{diag}{\left({\mathbf{\gamma }}_{{\Delta }^D}\right)}^{-1}({{\mathsf{\Sigma }}_X}_{{\Delta }^D}-\frac{2-{\lambda }_{{\Delta }^D}}{2{\theta }_{{\Delta }^D}}{\mathbf{\nu }}_{{\Delta }^D}^{\top }{\mathbf{\nu }}_{{\Delta }^D})\mathrm{diag}{\left({\mathbf{\gamma }}_{{\Delta }^D}\right)}^{-1}$, find the closest positive definite matrix and apply the denoising technique to estimate ${\mathsf{\Sigma }}_{{\Delta }^D}$.

3.3. Simulation

To conduct simulation with the calibrated model, we follow the procedures specified below.

Step 1: Simulate S sample paths of standard deviation of N assets for T days with the fitted model.

Step 2: Simulate S sample paths of a Markov chain ${\Delta }^D$ with the market’s transition matrix. Calculate the total number of each regime in all sample paths.

Step 3: Draw i.i.d. samples from stdMNTS distribution in each regime. The number of samples is determined by the calculation in Step 2. To draw one sample, we follow the procedure below.

Step 3.1: Sample from multivariate normal distribution N(0,${\mathsf{\Sigma }}_{{\Delta }^D}$).

Step 3.2: Sample from the tempered stable subordinator T with parameters ${\lambda }_{{\Delta }^D}$ and ${\theta }_{{\Delta }^D}$.

Step 3.3: Calculate $\mathbf{X}={\mathbf{\nu }}_{{\Delta }^D}(T-1)+\sqrt{T}({\mathbf{\gamma }}_{{\Delta }^D}\circ {\mathbf{\xi }}_{{\Delta }^D})$.

Step 4: Multiply the standard deviation in Step 1 with standard joint innovations in Step 3 accordingly. Add regime-specific mean ${\mathbf{\eta }}_{{\Delta }^D}$ to get simulated asset returns.

4. Portfolio Optimization

In this section, we discuss issues on portfolio optimization.

4.1. Problem Formulation

Conditioned on given sample path(s), classical portfolio optimization is formulated as a trade-off between risk and return

$$\begin{array}{cc}\hfill \underset{x}{min}& W\left(\mathcal{R}\right(x,t,\omega \left)\right|t\in [0,T],\omega \in \mathsf{\Omega })\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& E\left[\mathcal{R}(x,T,\omega )|\omega \in \mathsf{\Omega }\right]\ge d,\hfill \\ \hfill & x\in V,\hfill \end{array}$$

where d is the benchmark return, $\mathcal{R}(x,t,\omega )$ is the portfolio return at time $t\in [0,T]$ in given sample paths $\omega \in \mathsf{\Omega }$, $W\left(\mathcal{R}\right(x,t,\omega \left)\right|t\in [0,T])$ is a risk measure of $\mathcal{R}(x,t,\omega )$ for time $t\in [0,T]$ estimated based on given sample paths $\mathsf{\Omega }$, and $x=(x_1,\dots ,x_n)$ is the allocation vector of n assets. A typical setting of V is $\left\{x\right|{\sum }_{i=1}^N{x}_i=1,x_i\ge 0\}$, meaning that short selling is not allowed.

4.2. Discussion on Risk Measures and Constraint

If the historical sample paths of asset returns are used as input to the CVaR optimization, we minimize the CVaR of daily portfolio return from time $t=1$ to $t=T$. In the setting where we use simulation as input, we can choose to minimize the CVaR of daily returns from time $t=1$ to $t=T$ or CVaR of compounded returns at time $t=T$. We can set constraints on average daily returns from time $t=1$ to $t=T$ or compounded portfolio return at time $t=T$. Notice that the daily return in the simulation of GARCH model is not i.i.d. Furthermore, it is reasonable to focus on the the return at the final time. Thus, we choose to minimize the CVaR of return at time $t=T$ instead of return from time $t=1$ to $t=T$ subject to constraint on compounded portfolio return at time $t=T$.

The definition of drawdown needs some clarification as well. The portfolio drawdown at time t is defined in (13).

$${\mathrm{DD}}_{m,s}:=\underset{j\in \{0,1,2,\cdots ,m\}}{max}Q(x,t_j,{\omega }_s)-Q(x,t_m,{\omega }_s)=\sum _{i=h}^{m}r(x,t_i,w_s),$$

where $h=arg{max}_{j\in \{0,1,2,\cdots ,m\}}Q(x,t_j,{\omega }_s)$, is the sum of daily returns from its peak. It is somewhat unconventional at first sight compared with the absolute drawdown defined by

$${\mathrm{absDD}}_{m,s}:=\underset{j\in \{0,1,2,\cdots ,m\}}{max}P(x,t_j,{\omega }_s)-P(x,t_m,{\omega }_s)=P(x,t_h,{\omega }_s)\prod _{i=h}^m(1+r(x,t_i,w_s))$$

and the relative drawdown defined by

$${\mathrm{reDD}}_{m,s}:=\frac{{max}_{j\in \{0,1,2,\cdots ,m\}}P(x,t_j,{\omega }_s)-P(x,t_j,{\omega }_s)}{{max}_{j\in \{0,1,2,\cdots ,m\}}P(x,t_j,{\omega }_s)}=\prod _{i=h}^m(1+r(x,t_i,w_s)),$$

where $h=arg{max}_{j\in \{0,1,2,\cdots ,m\}}Q(x,t_j,{\omega }_s)$. CDaR is essentially an average of drawdowns that exceed the threshold. CDaR defined in (16) with ${\mathrm{DD}}_{m,s}$ satisfies good properties such as nonnegativity, insensitivity to constant shift, positive homogeneity and convexity (See Checkhlov et al. 2005), which gives great advantage in developing optimization techniques. CDaR can be defined in a similar way with ${\mathrm{absDD}}_{m,s}$ and ${\mathrm{reDD}}_{m,s}$.

It is possible to use the algorithm designed for CDaR defined by ${\mathrm{DD}}_{m,s}$ to minimize CDaR defined by ${\mathrm{absDD}}_{m,s}$. By substituting the input $r(x,t_i,w_s)$ with $P(x,t_i,{\omega }_s)-P(x,t_{i-1},{\omega }_s)$, the calculation of ${\mathrm{DD}}_{m,s}$ becomes

$$\sum _{i=h}^{m}P(x,t_i,{\omega }_s)-P(x,t_{i-1},{\omega }_s)={\mathrm{absDD}}_{m,s}$$

where $h=arg{max}_{j\in \{0,1,2,\cdots ,m\}}Q(x,t_j,{\omega }_s)$. However, it is obvious from the definitions that CDaR defined ${\mathrm{absDD}}_{m,s}$ and ${\mathrm{reDD}}_{m,s}$ do not satisfy basic axioms such as positive homogeneity.

Another drawback of absolute drawdown is that when we consider a long time period, drawdowns could vary greatly in absolute value but be close in relative value, the latter of which we are more concerned about. Uncompounded cumulative return is determined by rate of returns regardless of the price. This discussion justifies our definition of drawdown with ${\mathrm{DD}}_{m,s}$.

5. Empirical Study

In this section, we demonstrate the superiority of our approach with various in-sample and out-of-sample tests. We conduct model fitting, sample path simulation and portfolio optimization. In the in-sample test, we conduct a KS test on the marginal NTS distribution and marginal t-distribution to show that the NTS distribution has a much better fit on the innovations. We also provide the transition matrix and the correlation matrices of the joint innovations in each regime. In the out-of-sample test, we use various performance measures and visualizations for the rolling-window studies in recent years to show the outperformance of our portfolio.

For diversification purposes, we set the range of weights as $[0.01,0.15]$. A case study in Krokhmal et al. (2002) found that without constraints on weights, the optimal portfolio only consists of a few assets among hundreds. Investors often empirically categorize the market as bull, bear and sideways market. Haas et al. (2004) experiments with two and three regimes. We limit the number of regimes to be smaller than or equal to three. The Dip test is conducted on innovations to ensure unimodal innovation distribution. Starting fitting with three regimes, a model with fewer regimes is used instead when the p-value is lower than $0.1$ or have a regime lasting shorter than 40 days in total. Preliminary tests show that without these selection criteria, the model would have a risk of overfit in some cases. For example, the innovation could be multimodal, and the regime switches hundreds of times within the tested period. The phenomena is even more pronounced when the number of regimes is higher than three. To be concise in notation, we denote the portfolios in a convenient manner. For example, 0.9-CDaR portfolio denotes the optimal portfolio derived from the optimization with 0.9-CDaR as the risk measure. The standard deviation optimal portfolio is the optimal portfolio in the classical Markowitz model that maximizes Sharpe ratio, which we also denoted by mean-variance (MV) optimal portfolio. The word optimal is sometimes omitted. Since we simulate returns for 10 days, all CVaR with varying confidence levels are calculated with the simulated cumulative returns on the 10-th day. We omit the time and denote the optimal portfolio simply as CVaR portfolio.

5.1. Data

The data are comprised of the adjusted daily price of DJIA index, 29 of the constituents of DJIA Index and three ETFs of bond, gold and shorting index (Tickers: TMF, SDOW and UGL) from January 2010 to September 2020. One constituent DOW is removed because it is not included in the index throughout the time period. Since we use 1764 trading days’ (about 7 years) data to fit the model, the time period for out-of-sample rolling-window test is from January 2017 to September 2020, which includes some of the most volatile periods in history caused by the a pandemic and some trade tensions. By reweighting the index constituents, it can be regarded as an enhanced index portfolio.

5.2. Performance Measures

Various performance measures are discussed in Biglova et al. (2004), Rachev et al. (2008), Cheridito and Kromer (2013). Similar to Markowitz Model, an efficient frontier can be derived by changing the value of constraint. Each portfolio on the efficient frontier is optimal in the sense that no portfolio has a higher return with the same value of risk measure. In this paper, we refer to the portfolio with the highest performance ratio $\frac{E\left[R(x,T,·)\right]}{W\left(R\right(x,t,·\left)\right|t\in [0,T])}$ as an optimal portfolio and the others on the efficient frontier as suboptimal portfolios. Note that $E\left[R(x,T,·)\right]$ is ${\mathrm{CVaR}}_0\left(R(x,T,·)\right)$. Thus, the mean-CVaR ratio can be regarded as a special case of the Rachev ratio $\frac{{\mathrm{CVaR}}_{{\eta }_1}\left(R(x,T,·)\right)}{{\mathrm{CVaR}}_{{\eta }_2}\left(R(x,T,·)\right)}$. We test with ${\eta }_1={\eta }_2=0.1$ to focus on the distribution tails.

5.3. In-Sample Tests

We present various in-sample test results. The time period during which 2-regime model is first selected is used in fitting is from 30 April 2010 to 5 May 2017. 100 days is classified as regime 1, and 1665 days is classified as regime 2. The standard deviation of the DJIA index is $0.00217$ in regime 1 and $0.00801$ in regime 2. Regime 2 is more volatile than regime 1. One might expect that the market is in a less volatile regime longer and has a spiked standard deviation during a short but much more volatile regime. It is not the case in this tested period according to our model.

5.3.1. Kolmogorov–Smirnov Tests on Marginal Distributions

We conduct a Kolmogorov–Smirnov test (KS test) on the marginal distributions of the joint innovation of MRS-MNTS-GARCH model in the tested period. The results are presented in

Table 1. Kolmogorov–Smirnov tests on marginal distributions of joint innovation.

	Regime	Beta	KS Statistics of NTS	p-Value	KS Statistics of t	p-Value
AAPL	1	2.048	0.129	0.067	0.125	0.081
	2	0.038	0.018	0.624	0.055	0
AMGN	1	1.133	0.054	0.917	0.058	0.875
	2	0.139	0.012	0.955	0.045	0.002
AXP	1	−0.513	0.022	1	0.027	1
	2	−0.021	0.022	0.392	0.06	0
BA	1	−1.129	0.016	1	0.039	0.996
	2	−0.075	0.009	0.999	0.041	0.007
CAT	1	1.681	0.083	0.476	0.09	0.369
	2	0.069	0.012	0.973	0.043	0.004
CRM	1	−0.861	0.041	0.993	0.045	0.982
	2	0.095	0.032	0.071	0.071	0
CSCO	1	1.297	0.084	0.454	0.066	0.75
	2	0.058	0.043	0.004	0.08	0
CVX	1	−1.523	0.064	0.781	0.063	0.791
	2	−0.055	0.019	0.563	0.029	0.11
DIS	1	−0.644	0.115	0.132	0.103	0.225
	2	0.025	0.014	0.878	0.053	0
GS	1	1.274	0.038	0.998	0.031	1
	2	−0.046	0.014	0.917	0.036	0.026
HD	1	−1.431	0.035	0.999	0.056	0.901
	2	0.025	0.01	0.995	0.047	0.001
HON	1	−1.496	0.064	0.783	0.064	0.776
	2	0.131	0.015	0.843	0.054	0
IBM	1	0.64	0.023	1	0.041	0.993
	2	−0.019	0.019	0.587	0.057	0
INTC	1	−1.551	0.045	0.981	0.085	0.443
	2	0.061	0.008	1	0.044	0.003
JNJ	1	1.87	0.055	0.908	0.109	0.174
	2	−0.013	0.008	1	0.041	0.007
JPM	1	1.857	0.023	1	0.083	0.464
	2	0.026	0.008	1	0.041	0.007
KO	1	−1.674	0.072	0.656	0.078	0.546
	2	0.013	0.009	0.999	0.047	0.001
MCD	1	−1.956	0.067	0.742	0.092	0.343
	2	−0.005	0.022	0.396	0.061	0
MMM	1	0.658	0.019	1	0.025	1
	2	−0.1	0.014	0.882	0.054	0
MRK	1	2.367	0.156	0.013	0.157	0.013
	2	0.076	0.013	0.928	0.046	0.002
MSFT	1	−2.175	0.176	0.003	0.165	0.007
	2	0.26	0.016	0.753	0.057	0
NKE	1	1.652	0.05	0.956	0.1	0.258
	2	0.079	0.022	0.416	0.059	0
PG	1	1.035	0.03	1	0.062	0.815
	2	0.12	0.015	0.824	0.055	0
TLT	1	−0.283	0.056	0.897	0.042	0.992
	2	−0.161	0.036	0.025	0.028	0.144
TMF	1	−0.273	0.055	0.907	0.041	0.993
	2	−0.152	0.035	0.036	0.028	0.135
TRV	1	−1.532	0.055	0.91	0.058	0.864
	2	−0.033	0.011	0.988	0.048	0.001
UGL	1	1.139	0.066	0.749	0.068	0.721
	2	−0.074	0.012	0.973	0.051	0
UNH	1	1.543	0.05	0.951	0.098	0.274
	2	0.054	0.012	0.962	0.051	0
V	1	1.45	0.044	0.985	0.05	0.955
	2	0.016	0.025	0.26	0.061	0
VZ	1	−0.406	0.115	0.131	0.109	0.17
	2	−0.043	0.017	0.716	0.033	0.052
WBA	1	1.011	0.026	1	0.048	0.966
	2	0.166	0.022	0.368	0.062	0
WMT	1	1.804	0.08	0.52	0.076	0.582
	2	−0.047	0.016	0.8	0.054	0

For NTS distribution, there are 5 p-values smaller than $0.05$. For t-distribution, there are 31 p-values smaller than $0.05$. Most of the rejections correspond to Regime 2, indicating that the marginal t-distribution is outperformed by marginal NTS distribution in describing the return in the more volatile regime. We can also observe that $\beta$ varies significantly in different regimes. Generally, for all assets, $\beta$ in regime 1 has larger absolute values.

. For comparison, we also fit multivariate t-distribution on the joint innovations and report the results of KS tests on the marginal distributions. The degree of freedom of multivariate $\mathit{t}$-distribution in regime 1 and regime 2 are $5.34$ and $5.68$, respectively. Note that the marginal distribution of stdMNTS distribution is still NTS distribution. Since we assume common parameters $\lambda$ and $\theta$ in MNTS distribution, we report the parameter $\beta$ of each asset in both regimes. The data are rounded to 3 decimals place.

5.3.2. Transition Matrix and Denoised Correlation Matrices of Joint Innovation

The transition matrix of the innovations is presented in

Table 2. Transition matrix of joint innovation.

	Regime 1	Regime 2
regime 1	0.8538	0.1462
regime 2	0.0301	0.9699

This table presents the transition matrix of the joint innovation of MRS-MNTS-GARCH model.

. The self-transition probability in two regimes are close, and both are high.

We also provide the denoised correlation matrices of stdMNTS innovations in two regimes in

Table 3. Denoised correlation matrices of joint innovation in two regimes.

	AAPL	AMGN	AXP	BA	CAT	CRM	CSCO	CVX	DIS	GS	HD	HON	IBM	INTC	JNJ	JPM	KO	MCD	MMM	MRK	MSFT	NKE	PG	SDOW	TMF	TRV	UGL	UNH	V	VZ	WBA	WMT
AAPL	1	$0.15$	$0.14$	$0.09$	$0.17$	$0.23$	$0.21$	$0.09$	$0.14$	$0.12$	$0.15$	$0.25$	$0.16$	$0.13$	$0.28$	$0.08$	$0.10$	$0.04$	$0.20$	$0.11$	$0.19$	$0.23$	$0.07$	$-0.31$	$-0.05$	$0.08$	$-0.01$	$0.04$	$0.34$	$0.22$	$0.21$	$0.21$
AMGN	$0.28$	1	$0.15$	$0.11$	$0.13$	$0.40$	$0.18$	$0.12$	$0.27$	$0.16$	$0.27$	$0.29$	$0.20$	$0.18$	$0.30$	$0.16$	$0.02$	$0.06$	$0.24$	$0.39$	$0.13$	$0.29$	$0.13$	$-0.33$	$-0.22$	$-0.05$	$-0.02$	$0.04$	$0.31$	$0.13$	$0.25$	$0.10$
AXP	$0.34$	$0.39$	1	$0.37$	$0.44$	$0.15$	$0.29$	$0.33$	$0.38$	$0.50$	$0.21$	$0.38$	$0.13$	$0.13$	$0.22$	$0.48$	$0.04$	$0.15$	$0.29$	$0.28$	$0.01$	$0.27$	$0.09$	$-0.55$	$-0.23$	$0.18$	$-0.16$	$0.30$	$0.16$	$0.15$	$0.21$	$0.16$
BA	$0.36$	$0.40$	$0.50$	1	$0.33$	$0.09$	$0.24$	$0.13$	$0.25$	$0.28$	$0.09$	$0.35$	$0.13$	$0.05$	$0.22$	$0.28$	$0.08$	$0.13$	$0.27$	$0.18$	$-0.01$	$0.20$	$0.08$	$-0.44$	$-0.14$	$0.13$	$-0.09$	$0.26$	$0.06$	$0.17$	$0.19$	$0.12$
CAT	$0.38$	$0.32$	$0.48$	$0.49$	1	$0.18$	$0.30$	$0.30$	$0.26$	$0.36$	$0.15$	$0.43$	$0.16$	$0.28$	$0.10$	$0.33$	$0.04$	$0.06$	$0.30$	$0.15$	$0.16$	$0.11$	$0.10$	$-0.49$	$-0.16$	$0.14$	$-0.04$	$0.17$	$0.11$	$0.23$	$0.10$	$0.21$
CRM	$0.31$	$0.30$	$0.37$	$0.34$	$0.36$	1	$0.23$	$0.14$	$0.36$	$0.18$	$0.32$	$0.29$	$0.29$	$0.24$	$0.23$	$0.20$	$0.09$	$0.17$	$0.30$	$0.22$	$0.28$	$0.32$	$0.07$	$-0.33$	$-0.18$	$0.003$	$-0.10$	$-0.14$	$0.39$	$0.15$	$0.23$	$-0.03$
CSCO	$0.35$	$0.28$	$0.37$	$0.39$	$0.40$	$0.38$	1	$0.27$	$0.36$	$0.27$	$-0.02$	$0.30$	$0.25$	$0.28$	$0.19$	$0.29$	$0.16$	$0.13$	$0.30$	$0.21$	$0.24$	$0.09$	$0.16$	$-0.42$	$-0.15$	$0.12$	$-0.09$	$0.10$	$0.22$	$0.19$	$0.19$	$0.24$
CVX	$0.32$	$0.34$	$0.46$	$0.46$	$0.58$	$0.31$	$0.38$	1	$0.22$	$0.30$	$0.10$	$0.27$	$0.17$	$0.22$	$0.22$	$0.32$	$0.14$	$0.08$	$0.34$	$0.13$	$0.11$	$0.16$	$0.17$	$-0.45$	$-0.30$	$0.14$	$-0.04$	$0.13$	$0.11$	$0.20$	$0.19$	$0.17$
DIS	$0.33$	$0.39$	$0.49$	$0.49$	$0.44$	$0.37$	$0.41$	$0.44$	1	$0.32$	$0.29$	$0.40$	$0.26$	$0.09$	$0.25$	$0.32$	$0.11$	$0.19$	$0.34$	$0.08$	$0.16$	$0.34$	$0.19$	$-0.48$	$-0.29$	$0.13$	$-0.17$	$0.11$	$0.34$	$0.18$	$0.21$	$0.10$
GS	$0.34$	$0.39$	$0.55$	$0.48$	$0.50$	$0.35$	$0.41$	$0.47$	$0.49$	1	$0.16$	$0.31$	$0.09$	$0.21$	$0.10$	$0.64$	$-0.08$	$0.01$	$0.26$	$0.12$	$0.09$	$0.25$	$-0.01$	$-0.49$	$-0.34$	$0.10$	$-0.30$	$0.30$	$0.20$	$0.004$	$0.16$	$-0.04$
HD	$0.34$	$0.39$	$0.48$	$0.44$	$0.41$	$0.34$	$0.36$	$0.41$	$0.48$	$0.43$	1	$0.30$	$0.27$	$0.12$	$0.12$	$0.24$	$-0.04$	$0.14$	$0.29$	$0.13$	$0.18$	$0.29$	$0.04$	$-0.38$	$-0.23$	$0.13$	$-0.17$	$0.05$	$0.29$	$0.21$	$0.17$	$0.17$
HON	$0.41$	$0.43$	$0.56$	$0.63$	$0.64$	$0.42$	$0.46$	$0.57$	$0.57$	$0.55$	$0.54$	1	$0.27$	$0.22$	$0.28$	$0.37$	$0.21$	$0.24$	$0.47$	$0.17$	$0.27$	$0.28$	$0.31$	$-0.60$	$-0.31$	$0.31$	$-0.18$	$0.16$	$0.33$	$0.21$	$0.36$	$0.18$
IBM	$0.36$	$0.34$	$0.43$	$0.46$	$0.46$	$0.35$	$0.44$	$0.47$	$0.43$	$0.44$	$0.41$	$0.52$	1	$0.28$	$0.03$	$0.11$	$0.08$	$0.01$	$0.30$	$0.12$	$0.30$	$0.14$	$0.14$	$-0.39$	$-0.13$	$0.04$	$-0.10$	$-0.07$	$0.31$	$0.22$	$0.14$	$0.13$
INTC	$0.37$	$0.38$	$0.44$	$0.44$	$0.46$	$0.37$	$0.45$	$0.45$	$0.43$	$0.43$	$0.41$	$0.50$	$0.48$	1	$0.06$	$0.24$	$0.10$	$-0.04$	$0.27$	$0.17$	$0.45$	$0.08$	$0.16$	$-0.31$	$-0.01$	$-0.04$	$0.01$	$-0.02$	$0.22$	$0.21$	$0.12$	$0.09$
JNJ	$0.28$	$0.52$	$0.45$	$0.48$	$0.41$	$0.27$	$0.37$	$0.49$	$0.45$	$0.44$	$0.44$	$0.53$	$0.45$	$0.41$	1	$0.09$	$0.18$	$0.17$	$0.33$	$0.36$	$0.04$	$0.21$	$0.25$	$-0.37$	$-0.06$	$0.02$	$0.04$	$0.10$	$0.34$	$0.16$	$0.21$	$0.17$
JPM	$0.35$	$0.40$	$0.59$	$0.51$	$0.52$	$0.37$	$0.43$	$0.52$	$0.52$	$0.77$	$0.46$	$0.59$	$0.47$	$0.46$	$0.49$	1	$-0.04$	$0.09$	$0.30$	$0.15$	$0.14$	$0.24$	$0.05$	$-0.51$	$-0.35$	$0.24$	$-0.31$	$0.24$	$0.20$	$0.10$	$0.18$	$0.08$
KO	$0.27$	$0.35$	$0.39$	$0.42$	$0.35$	$0.28$	$0.33$	$0.46$	$0.43$	$0.33$	$0.40$	$0.47$	$0.43$	$0.40$	$0.48$	$0.40$	1	$0.23$	$0.10$	$0.13$	$0.21$	$0.06$	$0.35$	$-0.23$	$-0.02$	$0.40$	$0.06$	$0.05$	$0.17$	$0.26$	$0.31$	$0.14$
MCD	$0.30$	$0.35$	$0.35$	$0.39$	$0.34$	$0.33$	$0.33$	$0.38$	$0.38$	$0.33$	$0.42$	$0.42$	$0.40$	$0.35$	$0.44$	$0.35$	$0.47$	1	$0.26$	$0.09$	$0.02$	$0.15$	$0.09$	$-0.25$	$-0.16$	$0.35$	$-0.14$	$0.14$	$-0.01$	$0.02$	$0.24$	$0.07$
MMM	$0.39$	$0.42$	$0.52$	$0.56$	$0.57$	$0.37$	$0.45$	$0.55$	$0.53$	$0.51$	$0.50$	$0.67$	$0.51$	$0.50$	$0.56$	$0.55$	$0.48$	$0.45$	1	$0.17$	$0.13$	$0.31$	$0.21$	$-0.56$	$-0.30$	$0.29$	$-0.20$	$0.11$	$0.27$	$0.20$	$0.23$	$0.21$
MRK	$0.24$	$0.47$	$0.40$	$0.37$	$0.35$	$0.25$	$0.32$	$0.41$	$0.38$	$0.39$	$0.35$	$0.43$	$0.37$	$0.37$	$0.52$	$0.43$	$0.37$	$0.36$	$0.44$	1	$0.05$	$0.06$	$0.18$	$-0.36$	$-0.06$	$-0.03$	$0.05$	$0.27$	$0.24$	$0.16$	$0.18$	$0.22$
MSFT	$0.37$	$0.35$	$0.43$	$0.40$	$0.43$	$0.37$	$0.46$	$0.42$	$0.42$	$0.43$	$0.40$	$0.49$	$0.48$	$0.53$	$0.40$	$0.45$	$0.39$	$0.39$	$0.47$	$0.34$	1	$0.14$	$0.17$	$-0.28$	$0.06$	$0.03$	$0.0002$	$-0.04$	$0.31$	$0.29$	$0.21$	$0.12$
NKE	$0.30$	$0.32$	$0.41$	$0.41$	$0.36$	$0.38$	$0.33$	$0.33$	$0.41$	$0.36$	$0.45$	$0.48$	$0.35$	$0.35$	$0.36$	$0.39$	$0.37$	$0.43$	$0.43$	$0.31$	$0.37$	1	$0.12$	$-0.36$	$-0.18$	$0.12$	$-0.13$	$-0.06$	$0.33$	$0.07$	$0.21$	$0.14$
PG	$0.27$	$0.39$	$0.36$	$0.38$	$0.31$	$0.24$	$0.33$	$0.42$	$0.40$	$0.33$	$0.38$	$0.42$	$0.40$	$0.36$	$0.52$	$0.38$	$0.54$	$0.40$	$0.47$	$0.40$	$0.38$	$0.32$	1	$-0.32$	$0.01$	$0.15$	$0.01$	$0.08$	$0.31$	$0.19$	$0.33$	$0.21$
SDOW	$-0.48$	$-0.54$	$-0.67$	$-0.70$	$-0.70$	$-0.47$	$-0.56$	$-0.71$	$-0.68$	$-0.68$	$-0.64$	$-0.79$	$-0.69$	$-0.62$	$-0.68$	$-0.73$	$-0.61$	$-0.57$	$-0.78$	$-0.57$	$-0.61$	$-0.54$	$-0.58$	1	$0.35$	$-0.33$	$0.20$	$-0.34$	$-0.42$	$-0.34$	$-0.35$	$-0.30$
TMF	$-0.27$	$-0.24$	$-0.38$	$-0.33$	$-0.37$	$-0.27$	$-0.28$	$-0.34$	$-0.36$	$-0.42$	$-0.31$	$-0.38$	$-0.32$	$-0.30$	$-0.32$	$-0.45$	$-0.25$	$-0.24$	$-0.37$	$-0.26$	$-0.28$	$-0.29$	$-0.23$	$0.45$	1	$-0.27$	$0.41$	$-0.10$	$-0.18$	$-0.01$	$-0.10$	$0.03$
TRV	$0.30$	$0.39$	$0.46$	$0.46$	$0.45$	$0.31$	$0.38$	$0.48$	$0.49$	$0.51$	$0.45$	$0.54$	$0.45$	$0.42$	$0.50$	$0.56$	$0.48$	$0.42$	$0.56$	$0.39$	$0.41$	$0.38$	$0.47$	$-0.68$	$-0.30$	1	$-0.19$	$0.16$	$0.03$	$0.11$	$0.23$	$0.15$
UGL	$0.05$	$-0.06$	$-0.01$	$-0.01$	$0.10$	$-0.03$	$-0.01$	$0.14$	$-0.05$	$-0.05$	$-0.06$	$0.01$	$0.001$	$-0.01$	$-0.001$	$-0.07$	$0.01$	$-0.01$	$0.02$	$-0.02$	$-0.02$	$-0.07$	$0.03$	$-0.02$	$0.16$	$-0.01$	1	$-0.05$	$-0.22$	$0.09$	$-0.04$	$0.04$
UNH	$0.30$	$0.40$	$0.37$	$0.38$	$0.33$	$0.26$	$0.30$	$0.38$	$0.38$	$0.39$	$0.38$	$0.44$	$0.38$	$0.33$	$0.42$	$0.42$	$0.34$	$0.31$	$0.41$	$0.36$	$0.35$	$0.33$	$0.30$	$-0.54$	$-0.25$	$0.41$	$-0.002$	1	$-0.03$	$0.09$	$0.13$	$0.07$
V	$0.32$	$0.37$	$0.51$	$0.43$	$0.38$	$0.36$	$0.34$	$0.39$	$0.43$	$0.44$	$0.41$	$0.49$	$0.40$	$0.37$	$0.42$	$0.46$	$0.35$	$0.36$	$0.46$	$0.35$	$0.38$	$0.39$	$0.34$	$-0.60$	$-0.31$	$0.43$	$-0.02$	$0.34$	1	$0.16$	$0.17$	$0.13$
VZ	$0.26$	$0.36$	$0.39$	$0.35$	$0.36$	$0.23$	$0.30$	$0.42$	$0.40$	$0.34$	$0.37$	$0.44$	$0.38$	$0.36$	$0.47$	$0.40$	$0.46$	$0.37$	$0.46$	$0.38$	$0.36$	$0.31$	$0.43$	$-0.56$	$-0.21$	$0.46$	$0.01$	$0.26$	$0.33$	1	$0.18$	$0.27$
WBA	$0.25$	$0.36$	$0.33$	$0.36$	$0.31$	$0.28$	$0.26$	$0.32$	$0.34$	$0.33$	$0.35$	$0.41$	$0.33$	$0.31$	$0.41$	$0.35$	$0.33$	$0.31$	$0.38$	$0.34$	$0.33$	$0.29$	$0.36$	$-0.48$	$-0.22$	$0.36$	$-0.01$	$0.32$	$0.29$	$0.30$	1	$0.19$
WMT	$0.20$	$0.30$	$0.30$	$0.32$	$0.24$	$0.18$	$0.28$	$0.28$	$0.33$	$0.28$	$0.39$	$0.37$	$0.30$	$0.27$	$0.40$	$0.32$	$0.38$	$0.34$	$0.36$	$0.32$	$0.31$	$0.30$	$0.41$	$-0.47$	$-0.17$	$0.38$	$-0.02$	$0.29$	$0.27$	$0.37$	$0.30$	1

This table presents the denoised correlation matrices of joint innovations in two regimes. The upper triangle corresponds the regime 1 while the lower triangle corresponds to regime 2 which is more volatile. The correlations in regime 2 are generally higher than those in regime 1.

. Since there are two regimes identified in this period, we combine the two correlation matrices for easier comparison. The upper triangle of the matrix is the upper triangle of correlation matrix in regime 1, while the lower triangle of the matrix is the lower triangle of correlation matrix in regime 2. We find that the innovations have distinctively different correlation matrices in two regimes, validating the regime-switching assumption. The innovations are highly correlated in regime 2, which corresponds to more volatile periods.

5.4. Out-of-Sample Tests

We conduct out-of-sample tests to demonstrate the outperformance of the optimal portfolios. We use 3 types of risk measures with different confidence levels in the optimization, namely, maximum drawdown, 0.7-CDaR, 0.3-CDaR, average drawdown, 0.5-CVaR, 0.7-CVaR, 0.9-CVaR and standard deviation. Note that 0-CVaR is equal to the expected return, which does not make sense as a risk measure. Minimizing $\alpha$-CVaR with small $\alpha$ means that we include the positive part of the return distribution as risk. This is a drawback of standard deviation, and thus we set $\alpha$ at 0.5, 0.7, and 0.9 to demonstrate the superiority.

The rolling window technique is employed with 10-day forward-moving time window. In each time window, we simulate 10,000 sample paths of length 10 with the fitted model. The simulation is used as input to portfolio optimization. The portfolio is held for 10 trading days before rebalancing. The portfolio optimization is performed with software Portfolio Safeguard (PSG) with precoded numerical optimization procedures.

5.4.1. Optimal Portfolios

We report the performance of the optimal portfolios in various forms. Due to a large number of risk measures, we sometimes only present representative ones in graphs for better visualization. The other optimal portfolios have similar results.

Performance Measures

In

Table 4. Performance ratios of optimal portfolios with different risk measures.

	$\frac{\mathbf{Mean}\mathbf{Return}}{0\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.3\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.7\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{1\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.5\text{-}\mathbf{CVaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.7\text{-}\mathbf{CVaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.9\text{-}\mathbf{CVaR}}$	$\mathbf{Sharpe}\mathbf{Ratio}$	$\mathbf{Rachev}\mathbf{Ratio}$
0-CDaR	$0.072$	$0.051$	$0.027$	$0.008$	$0.258$	$0.156$	$0.077$	$0.117$	$0.965$
0.3-CDaR	$0.074$	$0.052$	$0.027$	$0.008$	$0.262$	$0.159$	$0.078$	$0.119$	$0.968$
0.7-CDaR	$0.073$	$0.052$	$0.027$	$0.008$	$0.260$	$0.158$	$0.077$	$0.118$	$0.965$
1-CDaR	$0.075$	$0.053$	$0.028$	$0.007$	$0.259$	$0.157$	$0.078$	$0.117$	$1.012$
0.5-CVaR	$0.076$	$0.054$	$0.028$	$0.007$	$0.259$	$0.158$	$0.076$	$0.120$	$1.002$
0.7-CVaR	$0.072$	$0.051$	$0.027$	$0.007$	$0.257$	$0.157$	$0.077$	$0.121$	$1.000$
0.9-CVaR	$0.071$	$0.050$	$0.027$	$0.007$	$0.243$	$0.147$	$0.072$	$0.116$	$0.986$
Standard Deviation	$0.064$	$0.045$	$0.024$	$0.006$	$0.238$	$0.145$	$0.071$	$0.109$	$0.966$
DJIA	$0.011$	$0.008$	$0.004$	$0.001$	$0.064$	$0.040$	$0.019$	$0.033$	$0.855$
Equal Weight	$0.032$	$0.023$	$0.011$	$0.003$	$0.136$	$0.083$	$0.037$	$0.065$	$0.900$

The column names are different ratios used to measure the performance. The row names show the risk measures used in the portfolio optimization. The DJIA index and the equally weighed portfolio are included as benchmarks.

, we measure the performance of the optimal portfolios by various performance ratios. The row names are the risk measures used in portfolio optimization. The columns are the performance measures. For example, the 0.3-CDaR portfolio is optimized every 10 days to maximize the ratio of expected cumulative return over 0.3-CDaR. Note that 0-CDaR is the average drawdown while 1-CDaR is the maximum drawdown. The risk measure standard deviation is studied in the classical Markowitz model, of which the optimization is also called MV optimization.

We can observe that the optimal portfolios of CDaR and CVaR measures have very close performance regardless of performance measures. The confidence level has a marginal impact on the ratios. Considering all the ratios, 0.3-CDaR portfolio and 0.5-CVaR portfolio are the best among the ten portfolios. 0.3-CDaR portfolio outperforms the others in all return-CVaR ratios, while 0.5-CVaR portfolio outperforms the others in return-CDaR ratios, slightly but consistently. All of the optimal portfolios of CDaR and CVaR measures consistently outperform the MV optimal portfolio, which significantly outperforms the DJIA and the equally weighted portfolio.

Log Cumulative Return

In Figure 1, we plot the time series of log cumulative returns of the optimal portfolios with different risk measures in out-of-sample tests to visualize the performance. The legend shows colors that represent different risk measures used in each portfolio optimization. We use log base 10 compounded cumulative returns as a vertical axis so that the scale of relative drawdown can be easily compared in graphs by counting the grids. The performance with different risk measures is reported separately. The labels in the legend show the confidence level of the risk measure. The names of the subfigures indicate the risk measures used in the optimization. DJIA index and the equally weighted portfolio are included in each graph for comparison.

As can be observed in the graphs, all the optimal portfolios follow a similar trend, though they differ in overall performance. They tend to alleviate both left tail and right tail events. CVaR optimal portfolios have higher cumulative returns than CDaR optimal portfolios. The confidence level has only a marginal impact. CVaR optimal portfolios with confidence levels 0.5 and 0.7 are almost indistinguishable. CDaR optimal portfolios with confidence levels 0.5, 0.7 and 0.9 are almost indistinguishable.

Relative Drawdown

Since the CVaR and CDaR risk measures only concern the tail behavior of a portfolio, we plot the relative drawdown of the optimal portfolios in out-of-sample testing period in Figure 2 to demonstrate their ability to alleviate extreme left tail events. Due to a large number of optimal portfolios with different risk measures, we only plot 0-CDaR, 0.5-CVaR and MV optimal portfolios here for better visualization. The others have similar results. The 0-CDaR, 0.5-CVaR and MV optimal portfolios have significantly smaller relative drawdown most of the time in the out-of-sample test, especially in year 2020. 0-CDaR and 0.5-CVaR optimal portfolios are slightly better than the MV optimal portfolio.

Portfolio Return Distribution

In Figure 3, we plot the kernel density estimation on the return distribution. To better compare the tails of the return distribution, we also plot the log base 10 density. We can easily observe that the optimal portfolios are very close and alleviate both right and left tail events. Unlike standard deviation (MV) optimal portfolio, 0.3-CDaR and 0.5-CVaR portfolio have thinner tails on both sides than DJIA and the equally weighted portfolio. 0.3-CDaR and standard deviation portfolio have a higher peak than the equally weighted portfolio. Among the three optimal portfolios, the 0.5-CVaR portfolio has the thinnest left tail and a fatter right tail than the other two optimal portfolios.

Allocation

In Figure 4, we report the time series of optimal weights of the 0-CDaR optimal portfolio as a representative case. The DJIA constituents are aggregated for better visualization. It is reasonable that when the total weight on DJIA constituents is high, the weights on shortselling ETF is low. The weights on ETFs hit the lower and upper bounds 0.01 and 0.15 in many periods, indicating that the constraints on weights are active in the optimizations.

5.4.2. Robustness to Suboptimality in Optimization

So far, we have been studying the optimal portfolios that maximize performance ratios. In this section, we study the sensitivity of the performance by comparing suboptimal portfolios with the optimal ones. We show that the portfolio with tail risk measure is more robust to suboptimality in optimization than the mean-variance portfolio. The performance ratios of CDaR and CVaR optimal portfolios are shown to be not sensitive to changing the value of return constraint, while the performance ratios of mean-variance optimal portfolio shows significant change. The graphs of log cumulative return show that the suboptimal portfolios follow a similar trend, sometimes overlapping with each other.

We use 9 portfolios with different constraints on return to approximate the efficient frontier for each risk measure. For example, in each 0.3-CDaR optimization, we perform 10 optimizations with varying constraints on return. The portfolio with the highest return-risk ratio is the approximated optimal portfolio. We denote the portfolio that has n lower level of constraint on return as level Ln suboptimal portfolio, the one that has n higher level of constraint on return as level Hn suboptimal portfolio. When a certain optimization is unfeasible, we set the allocation the same as that of a lower level. The constraints on return are (0.002, 0.010, 0.020, 0.030, 0.035, 0.040, 0.045, 0.050, 0.060, 0.065). For example, when the portfolio with constraint 0.050 on return has the highest performance ratio among the 10 portfolios, it is labeled as the optimal portfolio. Accordingly, the portfolio with one-level higher constraint 0.060 is labeled H1. In this case, H2 to H4 have the same constraint 0.065, since there is no optimization performed with higher constraint. Note that at each rebalance, the optimal portfolio changes, so do portfolios of other levels, i.e., constraint 0.060 may no longer be optimal. The performance of each label is measured with portfolios of corresponding level determined at each rebalance.

We report suboptimal portfolios with risk measures of 0-CDaR, 0.5-CVaR and standard deviation in Figure 5 and

Table 5. Performance ratios of suboptimal portfolios with different risk measures.

	$\frac{\mathbf{Mean}\mathbf{Return}}{0\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.3\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.7\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{1\text{-}\mathbf{CDaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.5\text{-}\mathbf{CVaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.7\text{-}\mathbf{CVaR}}$	$\frac{\mathbf{Mean}\mathbf{Return}}{0.9\text{-}\mathbf{CVaR}}$	$\mathbf{Sharpe}\mathbf{Ratio}$	$\mathbf{Rachev}\mathbf{Ratio}$
0-CDaR
L4	$0.061$	$0.043$	$0.023$	$0.008$	$0.239$	$0.146$	$0.074$	$0.113$	$0.965$
L3	$0.058$	$0.041$	$0.021$	$0.007$	$0.232$	$0.142$	$0.071$	$0.108$	$0.959$
L2	$0.061$	$0.043$	$0.023$	$0.007$	$0.236$	$0.144$	$0.072$	$0.110$	$0.960$
L1	$0.062$	$0.044$	$0.023$	$0.007$	$0.236$	$0.144$	$0.072$	$0.110$	$0.959$
optimal	$0.072$	$0.051$	$0.027$	$0.008$	$0.258$	$0.156$	$0.077$	$0.117$	$0.965$
H1	$0.074$	$0.052$	$0.028$	$0.008$	$0.251$	$0.151$	$0.073$	$0.115$	$0.968$
H2	$0.082$	$0.058$	$0.030$	$0.008$	$0.259$	$0.157$	$0.076$	$0.120$	$0.993$
H3	$0.087$	$0.061$	$0.032$	$0.009$	$0.267$	$0.161$	$0.078$	$0.124$	$0.995$
H4	$0.087$	$0.062$	$0.032$	$0.009$	$0.268$	$0.162$	$0.079$	$0.125$	$0.989$
0.5-CVaR
L4	$0.073$	$0.052$	$0.027$	$0.006$	$0.252$	$0.154$	$0.074$	$0.114$	$1.011$
L3	$0.074$	$0.052$	$0.027$	$0.007$	$0.253$	$0.155$	$0.075$	$0.115$	$1.011$
L2	$0.074$	$0.052$	$0.027$	$0.007$	$0.253$	$0.155$	$0.075$	$0.116$	$1.008$
L1	$0.075$	$0.053$	$0.027$	$0.007$	$0.256$	$0.156$	$0.076$	$0.118$	$1.006$
optimal	$0.076$	$0.054$	$0.028$	$0.007$	$0.259$	$0.158$	$0.076$	$0.120$	$1.002$
H1	$0.075$	$0.053$	$0.027$	$0.007$	$0.255$	$0.155$	$0.075$	$0.118$	$0.990$
H2	$0.075$	$0.053$	$0.027$	$0.007$	$0.256$	$0.156$	$0.075$	$0.119$	$0.987$
H3	$0.073$	$0.052$	$0.027$	$0.007$	$0.256$	$0.155$	$0.075$	$0.119$	$0.982$
H4	$0.071$	$0.050$	$0.026$	$0.007$	$0.254$	$0.154$	$0.074$	$0.117$	$0.977$
Standard Deviation
L4	$0.041$	$0.029$	$0.015$	$0.004$	$0.176$	$0.110$	$0.055$	$0.084$	$0.990$
L3	$0.042$	$0.030$	$0.016$	$0.004$	$0.177$	$0.111$	$0.056$	$0.083$	$1.005$
L2	$0.049$	$0.035$	$0.019$	$0.005$	$0.196$	$0.122$	$0.062$	$0.092$	$1.013$
L1	$0.059$	$0.042$	$0.022$	$0.006$	$0.225$	$0.138$	$0.069$	$0.103$	$1.005$
optimal	$0.064$	$0.045$	$0.024$	$0.006$	$0.238$	$0.145$	$0.071$	$0.109$	$0.966$
H1	$0.069$	$0.049$	$0.026$	$0.007$	$0.232$	$0.141$	$0.069$	$0.108$	$0.972$
H2	$0.071$	$0.050$	$0.027$	$0.007$	$0.232$	$0.141$	$0.069$	$0.110$	$0.983$
H3	$0.074$	$0.052$	$0.028$	$0.007$	$0.237$	$0.144$	$0.070$	$0.113$	$0.986$
H4	$0.073$	$0.052$	$0.027$	$0.007$	$0.234$	$0.142$	$0.069$	$0.112$	$0.980$
DJIA	$0.011$	$0.008$	$0.004$	$0.001$	$0.064$	$0.040$	$0.019$	$0.033$	$0.855$
Equal Weight	$0.032$	$0.023$	$0.011$	$0.003$	$0.136$	$0.083$	$0.037$	$0.065$	$0.900$

This table presents the performance ratios of suboptimal portfolios with different risk measures. The risk measures are shown in the row names: 0-CDaR, 0.5-CVaR and standard deviation. The constraint on return is adjusted from high to low (H4 to L4) to show the impact of suboptimality. Portfolio optimization with tail risk measures demonstrates smaller variation in out-of-sample performance ratios compared with mean-variance portfolio.

as representative cases. We find that the optimal portfolio is consistently too conservative for all risk measures and input. That is, with the same risk measure, suboptimal portfolios with a higher constraint on return achieve higher performance ratios than the optimal portfolio. From $H4$ to $L4$, the ratios do not show a pattern of first increasing the decreasing as expected. In other words, the out-of-sample efficient frontiers are usually non-convex with irregular shapes. Overall, the performance ratios of CDaR and CVaR optimal portfolios are shown to be not sensitive to changing the value of return constraint, while the performance ratios of mean-variance optimal portfolio shows significant change.

The portfolios have increasing returns as well as risk from L1 to H4. The realized paths follow a similar trend and rarely cross. The optimal portfolios, as expected, are in the medium part of all paths. The 0.5-CVaR suboptimal portfolios have almost identical performance. Some paths are not visible due to overlap. For example, the constraints on the return of H4 are sometimes not feasible, leading to the same performance with H3 portfolio.

6. Conclusions

We propose the MRS-MNTS-GARCH model that is sufficiently flexible to accommodate fat tails, skewness and regime switch in asset returns. The model is used to simulate sample paths of portfolio returns, which serves as input to portfolio optimization with tail risk measures. We conduct various in-sample and out-of-sample tests to demonstrate the effectiveness of our approach. In-sample tests show that the NTS distribution fits the innovations with high accuracy compared with t-distribution. Out-of-sample tests show that our approach significantly improves the performance of optimal portfolios measured by performance ratios and it successfully mitigates left tail risk. We also find that out-of-sample performance ratios of the portfolios with tail risk measures are more robust to suboptimality on the efficient frontier compared with mean-variance portfolio.