Empirical studies have found in the return of various financial instruments skewness and leptokurtoticasymmetry, and higher peak around the mean with fat Pareto-type tails. Normal distribution has long been recognized as insufficient to accommodate these stylized fact, relying on which could drastically underestimate the tail risk of a portfolio. However, few has adequately incorporate this fact into modeling and decision-making. The α-Stable distribution has been proposed to accommodate them. The lack of existence of moments of stable distributions could sometimes cause difficulties, to which tempered stable distribution presents a potential solution.
The class of tempered stable distributions is derived from the α-stable distribution by tempering tails. Recent findings from Kim et al. (2012) , Kurosaki and Kim (2018) , Anand et al. (2016) use the normal tempered stable(NTS) distribution and successfully model the stock returns with high accuracy. As is pointed out, the above mentioned properties not only exist in raw returns but also in GARCH-filtered residuals of time series model. This motivates us to use multivariate normal tempered stable distribution to accommodate the asymmetry, interdependence and fat tail of the joint innovations in our model. Other distributions like Pearson distribution and generalized hyperbolic distribution have also been proposed to this end, both having t distribution as a special case. The estimation on the degree of freedom parameter of t distribution in GARCH model is found to be unreliable. Study in Kim et al. (2011) finds that normal and t distribution are rejected and normal temped stable distribution is favored on the residuals of time series model. Another study in Shi and Feng (2016) found that tempered stable distribution is a better choice than t or GED distribution in MRS-GARCH innovation. While generalized hyperbolic distribution is very flexible, it has 5 parameters which might create difficulty in estimation. NTS distribution has only 3 parameters and is flexible enough to serve the purpose.
One of the drawbacks of GARCH model is the difficulty in dealing with volatility spikes, which could indicate that the market is switching within regimes. Various types of specification of Markov regime switching GARCH models has been proposed, many of which suggest that regime switching GARCH model achieves a better fit for empirical data. Naturally, the most general model allows all parameters to switch among regimes. However, as is shown in Henneke et al. (2011) , the sampling procedure in MCMC method for the estimation of such model is time-consuming and renders it improper for practical use. Recent developments in algorithm in Billioa et al. (2016) improves the estimation speed. In our model, we use the regime switching GARCH model specified in for simplicity. This model circumvents the path dependence problem in Markov model by specifying parallel GARCH models.
Its a recognized fact that the correlation of financial assets is time-varying. Specification of multivariate GARCH models with both regime-specific correlations and variance dynamics involves a balance between flexibility and tractability. The model in has been 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 Markov swiching multivariate NTS distribution.
Modern portfolio theory is usually framed as a trade-off between return and risk. Classical
Markowitz Model intends to find the portfolio with the highest Sharpe ratio. However, variance has been criticized for not containing enough information on the tail of distribution. Moreover, the correlation is not sufficient to describe the interdependence of asset return 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. However, there are several undesired properties that rendered it an insufficient criterion. First, its not a coherent risk measure due to lack of sub-additivity. Second, VaR is difficult to to optimize for non-normal distributions due to non-convex and non-smooth as a function of positions with multiple local extrema, which causes difficulty in developing efficient optimization techniques. Third, a single percentile is insufficient to describe the tail behavior of a distribution, which might lead to an underestimation of risk.
Theory and algorithm for portfolio optimization with CVaR measure is proposed in Krokhmal et al. (2001a) Rockafellar and Uryasev (2000) to address these issues. For continuous distributions, CVaR is defined as a conditional expectation of losses exceeding a certain VaR level. For discrete distributions, CVaR is defined as 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 measure can be found in Sarykalin et al. (2014) .
Drawdown has been a constant concern for investors. Its 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. Thus maximum drawdown is often used in evaluation of performance of a portfolio. As is mentioned in Checkhlov et al. (2005) , a client account of a Commodity Trading Advisor(CTA) will be issued a warning or shut down with a drawdown higher than 15% or longer than 2 years, even if small. Investors are highly unlikely to tolerate a drawdown greater than 50%. However, maximum drawdown only considers the worst case which may only occur under some very special circumstances. It's also very sensitive to the testing period and asset allocation. On the other hand, small drawdowns included in average drawdown are acceptable and might be caused by pure noise, of which minimization might not make sense. For instance, a Brownian motion would have drawdowns in a simulated sample path.
A relevant criterion to CVaR, conditional drawdown(CDaR) is proposed in Checkhlov et al. (2005) to address these concerns. While CVaR only concerns the distribution of return, CDaR concerns the sample path. CDaR is essentially a certain CVaR level of the drawdowns. By this, we overcome the drawbacks of average drawdown and maximum drawdown. CDaR not only take the depth of drawdowns into consideration, but also the length of them, corresponding to the concerns by investors we mentioned earlier. 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. Optimization method with constraint on CDaR has also been developed in relevant papers.
For an optimization procedure to lead to desired optimal result, the input is of critical importance. With careless input, the portfolio optimization could magnify the impact of error and lead to unreasonable result. Using historical sample path like in Krokhmal et al. (2001b) means that we assume what happened in the past will happen in the future, which is an assumption that needs careful examination. It's 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 sample paths. This motivates us to propose a model that incorporates fat tails to generate simulation as input for the optimization with tail risk measures.
To this end, we propose a Markov regime switching GARCH model with multivariate standard normal tempered stable innovations. In-sample tests are performed to examine the goodness of fit. Simulation is performed with the model to generate multiple sample paths. To demonstrate the effectiveness of the proposed model as well as the tail risk measures, we conduct an ou-ofsample study on the performance. The tested period includes some of the most volatile time in history caused by COVID-19 pandemic and international trade tensions.
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, methods for estimation and simulation. Section 4 is an empirical study on in-sample goodness of fit and out-of-sample performance in recent years.
The normal tempered stable (NTS) distribution is defined by the subordinated Gaussian distribution as follows. Let T be a strictly positive random variable defined by the characteristic function for λ ∈ (0, 2) and θ > 0
The n-dimensional NTS distributed random vector X = (X 1 , ..., X n ) is defined as
Where µ, ν ∈ R n , γ ∈ R n + , ξ ∼ N(0, Σ), γ • ξ = (γ 1 ξ 1 , ...., γ n ξ n ). T is the Tempered Stable subordinator with parameters λ and θ, which is independent of ξ. This is an infinite mixture of multivariate normal distribution. Multivariate NTS distribution specified above is denoted as MNTS(λ, θ, ν, γ, µ, Σ). A bold letter denotes a vector whenever used.
A standard multivariate NTS distribution has µ = (0, ...., 0),
2−λ , where ν i and γ i are the ith elements of ν and γ. This yields a distribution that has unit variance for conditional distribution. The covariance between X i and X j is given by
where diag(γ) is the diagonal matrix with γ i , i = 1, ..., n as diagonal elements. This formula is used for subsequent parameter estimation. The pdf of standard multivariate NTS distribution is given by
Setting the dimension n = 1 will give us univariate NTS distribution. The marginal distribution of an MNTS distribution is still an NTS distribution.
GARCH(p, q) model has been studied intensively as a model for volatility.
where r t is the return at time t, σ t is the variance at time t, η, ω, α i , β i are parameters. We will use GARCH(1,1) in the paper.
Before specifying our model, it would be clear if we first specify the univariate Markov regime switching model in .
where ∆ t is a Markov chain with finite state space S = 1, ..., k and an irreducible and primitive k × k transition matrix P, with element p ij = P (∆ t = j | ∆ t−1 = i), σ 2 t = (σ 2 1,t , ...., σ 2 k,t ), ω = (ω 1 , ...., ω k ), α = (α 1 , ...., α k ), β = (β 1 , ...., β k ), • denotes element-wise product. The stationary conditons require a definition of matrices In principle, regime-specific mean can be included in this model, though it may cause loss in property of zero auto-correlation.
Clarifying notations would facilitate the following illustration of rsisk measures. Denote the rate of return of asset n from time t − 1 to t of the ith sample path as r (i) tn , i = 1, ...., S, t = 1, ...., T, n = 1, ...., N . In this paper, return always refers to rate of return unless specified as cumulative return, R and r always denote return of some types, sample path refers to a time series of rate of return rather than asset price.
U (x) can be mean return or cumulative return of the portfolio at the end of the considered period, both of which are specified as linear constraints in the optimization problem. When we consider cumulative return, it can be either compounded or uncompounded. Denote the time series of cumulative portfolio return as
when uncompounded. Denote the probability of the ith sample path as
T . In our simulation study, all simulated sample paths have equal probability 1 S . Thus we don't explicitly show p (i) in subsequent discussion on multiple sample paths. (2000), we assume a continuous loss function f (x, y) where x is the asset allocation vector, y is the random variable of return, p(y) is the PDF of y. The probability of f (x, y) not exceeding a threshod ζ is
The α-VaR is defined as
The CVaR risk measure of level α, denoted as α−CVaR or CVaR α , is given by
It is also called average VaR due to the formula.
When we consider an accumulated portfolio return R P in the time series of a single sample path as f (x, y), y becomes a discrete distribution with point value {r 1P , ...., r T P } and p(y) has equal point probability.
To be clear about the computation of discrete case which we will consider, we follow Checkhlov et al. (2005) to introduce a function
where I is an indicator function that equals 1 when the condition in the bracket is true and 0 when it's false.
The inverse function of π R P (s) is defined as
CVaR in discrete case is thus given by
This is the formula we use when we have historical path as input. Since CVaR is not pathdependent, it may come to mind that the definition still applies In the case of multiple simulated sample paths R (i) P , i = 1, ...., N by binding all the simulated path into one. One of the drawback to calculate VaR and CVaR with historical data is the implicit assumption that the return are i.i.d. For the same reason, it's inappropriate to aggregate all the simulated returns in our case, since simulation with GARCH model clearly means that the simulated daily returns days ahead are not i.i.d. In the case of simulated sample paths, we should use only the T −day accumulated return of each sample path to calculate T −day VaR and CVaR, i.e. substitute R
As is explained above, it's not necessary to define CVaR for multiple sample paths in our case.
Nevertheless, the following definition facilitates subsequent definition of conditional drawdown.
A similar function, its inverse function and CVaR are respectively defined
Drawdown is a measurement of the price fall from the last peak. It can be defined in both absolute and relative way. CDaR α measures the 100 · α percent of largest drawdowns. The algorithm is developed for absolute drawdown since it has a close relation with CVaR and thus possesses some desired property in optimization. For each cumulative portfolio return series
T P ), we define its drawdown series
Conditional drawdown of level α, denoted as α−CDaR or CDaR α , is the α% largest percent of all drawdowns. Its definition is an application of CVaR in a dynamic case. We respectively define α-CDaR of single sample path and multiple sample paths simply by substituting with
This includes average drawdown and maximum drawdown as special cases where α equals 1 and 0.
To model a multivariate process, we assume the variance of individual asset evolves according to univariate regime-switching GARCH model independently with possibly different number of regimes, while the correlation is reflected separately in the joint standard residuals.
regime. In the model we use, however, the variance at time t is determined by the variance at time t − 1 within the regime. The joint distribution of innovations is modeled with a multivariate NTS distribution and evolves according to a Markov Chain. That is, the regimes of each asset and joint innovations evolve independently, avoiding the difficulty where all parameters switch regimes simultaneously. Intuitively, the correlation of assets change drastically when a market regime shift takes place. Thus here for simplicity we assume that the multivariate innovations follow the same path as the market which is denoted as ∆ M k .
Though flexible enough to accommodate the many stylized facts of asset return, the lack of an analytical form of NTS distribution presents some difficulties in estimation. Our estimation methodology is adapted from Kim et al. (2011) , Kim et al. (2012) , which is among the first works to incorporate NTS distribution in portfolio optimization.
We first fit univariate model on each asset and extract the residuals for subsequent estimation of MNTS distribution. The MNTS distribution has two tail parameters λ, θ and one skewness parameter vector ν. Common tail parameters λ, θ is assumed for individual constituents and estimated from the DJIA index in each regime respectively. This leaves the skewness parameters ν to be calibrated by MLE for each asset. We explicitly estimate the covariance matrix of the innovations in each regime to compute the skewness parameters. It's known that historical correlation matrix is unreliable. We find that the out-of-sample performance is enhanced with denoising techniques in Laloux et al. (1999) , Juliane and Korbinian (2005) .
The joint residuals are assumed to follow a MNTS Hidden Markov process. In our approach, the MNTS distribution is estimated within each market regime. Estimation of MNTS distribution with EM algorithm combined with FFT is studied in Bianchi et al. (2016) . It's probably feasible but time-consuming to fit the MNTS Hidden Markov Model on residuals with MCEM algorithm.
We find that univariate Markov switching model with more than 3 states often has one unstable state that lasts for a very short period and switches frequently and sometimes has one state with clearly bimodal residuals distribution. Thus it's desirable to limit the number of states smaller than 4, check unimodal distriution with Dip test and choose the one with highest BIC value.
To conduct simulation with the model, we first need to simulate the path of variance for each asset and the index with univariate Markov regime switching GARCH model. Then we simulate the multivariate innovations with MNTS distribution, whose regime switching probabilities follow the indexs. Finally, we multiple the standardized innovations with the variance and add the mean value to get simulated return.
With given sample path(s), classic portfolio optimization is formulated as a trade-off between risk and return min
where w denotes the given sample path(s), U (x) is portfolio return and M (x) is a risk measure,
x = (x 1 , ...., x n ) is the allocation vector of n assets. Since there are multiple sample paths in the case of simulation, the expected value is used. A typical setting of V is {x | N i=1 x i = 1, x i ≥ 0}, meaning that short selling is not allowed. Sharpe-like ratio can be defined by substituting standard deviation in Sharpe ratio with other risk measures.
Similar to Markowitz Model, efficient frontier can be derived by changing d. Each portfolio on the efficient frontier is optimal in the sense that no portfolio has higher return with the same value of risk measure. For clarification in this paper, we refer to the portfolio with the highest Sharpe-like ratio Ew(M (x|w)) Ew(U (x|w)) as optimal portfolio and the others on efficient frontier as suboptimal portfolios.
The algorithm developed in relevant papers aims to minimize absolute drawdown specified above. While the properties of absolute drawdown facilitates the computation, a natural concern to this approach is that when we consider a long period, drawdowns at different times may differ greatly in absolute value but close in relative value, of which the latter we may care more about.
Consider drawdown w
nP . Using uncompounded cumulative return in the setting could alleviate this problem, since the drawdown ξ k would be essentially sum of returns k j=n r et al. (2001b) found that without constraints on weight, the optimal portfolio only consists of a few assets among hundreds. We limit the number of states to be smaller than 4. Dip test is conducted on residuals to ensure unimodal residual distribution. When the p-value is lower than 0.1, model with fewer states is used instead. With this condition satisfied, the one with highest BIC is selected. To be concise in notation, we denote the portfolios in a convenient manner. e.g., 0.9-CDaR portfolio denotes the optimal portfolio derived from the optimization with 0.9-CDaR as risk measure.
The data comprises of the adjusted daily price of DJIA index, 29 of the constituents of DJIA Index and 3 ETFs TMF, SDOW and UGL, from January 2010 to September 2020. One constituent is removed for that it's not included in the index in all tested period. Since we use 1764 trading days' (about 7 years) data to fit the model, the actual tested period is from January 2017 to September 2020, which includes some of the most volatile periods in history caused by the COVID-19 pandemic and trade tensions.
We present an in-sample test result. The time period used in fitting is from 2013-07-19 to 2020-07-22. 2-state model is favored in this period. 1075 days is classified as state 1 and 690 days is classified as state 2.
To examine the goodness of fit of SMNTS distrbution on the joint residuals, we report in Table I the p value and KS statistics of KS test on residuals of each asset and marginal distribution of SMNTS distribution, which is still a standard NTS distribution. As a comparison, we fit multivariate Student's t distribution on the joint residuals and conducted KS test on marginal distribution similarly. We also report the β of each asset in different regimes.
The degree of freedom of multivariate Student's t distribution in 2 regimes is 5.6 and 16.1 respectively. The data are rounded to 3 decimals place.
We can observe that β varies significantly in different regimes. Generally, for all assets, β in regime 2 has much larger absolute values, indicating a significant skewness. NTS fits very well for all residuals in regime 1 due to its ability to model distribution with a high peak, skewness and fat tails, where t distribution are frequently rejected with p value close to 0. In regime 2, NTS distribution performs similarly to Student's t distribution.
2) Residuals: Transition matrix of the residuals is shown in Table II .
We provide the denoised correlation matrices of MNTS residuals in 2 regimes in Table III .
Since there are 2 regimes identified in this period , we stack the upper tridiagonal part of the matrix with elements in the 1st regime and the lower tridiagonal part with the 2nd together for easier comparison. We find that the residuals have distinctively different correlation matrices in 2 regimes, validating the regime switching assumption. The residuals are highly correlated in the 2nd regime that corresponds to more volatile periods with a lower self-transition probability.
Out-of-sample test with real market data is an informative way to test the effectiveness of a model. We use 3 types of risk measure with different confidence level 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 doesn't make sense as a risk measure. Minimizing α-CVaR with α smaller than 0.5 means that we include the right part of the return distribution with many positive values. This is a flaw of standard deviation and thus we set α at 0.5, 0.7,and 0.9 to demonstrate the superiority.
Rolling window technique is employed with 1764 days forward moving time window for biweekly estimation of the parameters. We simulate 1000 sample paths of length 10 every two weeks with the fitted model. The simulation is used as input for portfolio optimization. The portfolio is hold for two weeks until rebalance. The portfolio optimization is performed with software PSG with precoded numerical optimization procedure. The L A T E X code for tables is produced with Hlavac (2018).
We report the performance of the optimal portfolios in various forms.
a) Ratio Measurement: We use Sharpe-like ratios, i.e. the mean return of realized path divided by a risk measure of realized path, to measure performance of the optimal portfolios. The result is reported in Table IV . The row names denote the risk measure in the portfolio optimization. The columns are the performance measure.
The outperformance of our model and the two tail risk measures in tems of the ratios is apparent. The performance measures are consistent in the sense that a portfolio outperforming in one measure also outperforms in other measures as well in most cases. We find that 2 tail risk measures lead to similar performance. The confidence level only has a slight impact. Both CVaR and CDaR optimal portfolios slightly but consistently outperform standard deviation optimal portfolio.
. b) Accumulated Return: To visualize the performance, we report in Fig.1 As can be observed in the graphs, all the optimal portfolios follow a similar trend, though differ in overall performance. They tend to avoid both left tail and right tail events. For all CDaR optimal portfolios with different input, confidence level have only a small impact on the realized path, which often overlap. It's possibly because the simulation length is only 10 days and the drawdown distribution is far from continuous, leading to close calculation of drawdowns with different confidence level.
CDaR risk measures concern the tail behavior, we plot the relative drawdown series and the (1) Optimization with the simulated sample paths by our model leads to good optimal portfolios with higher performance than the index and equally weighted portfolio.
(2) Tail risk measure CDaR and CVaR is superior to standard deviation in optimization.
(3) Using our model and CDaR measure together leads to best performance, while model has a greater impact than the risk measure.
Due to the large number of optimal portfolios with different risk measures, we only study 0-CDaR, 0.5-CVaR and standard deviation optimal portfolios here. The others have similar results. As is shown in Fig.2 , the 0-CDaR, 0.5-CVaR and standard deviation optimal portfolios have significantly smaller relative drawdown in most period in the out-of-sample test, while 0-CDaR and 0.5-CVaR optimal portfolios are slightly better than standard deviation portfolio.
From the distribution in Fig.3 , we can see that the optimal portfolios are almost identical and successfully avoid both right and left tail events. They have thinner tails on both sides than DJIA and equally weighted portfolio and slightly skewed to the right. The optimal portfolios have lower peak than the equally weighted portfolio but still much higher peak than DJIA. We also report the series of allocation weights among the DJIA constituents and the ETFs in Fig.4 . It's reasonable that when the weights on DJIA constituents are high, the weights on shortselling ETF is low.
2) Suboptimal Portfolio: So far we have been studying the optimal portfolios that maximize Sharpe-like ratios. In this section, we study whether the optimal portfolio outperforms suboptimal ones.
We use 9 portfolios with different constraint on return to approximate the efficient frontier for each risk measure. For example, in each 0.3-CDaR optimization, we perform 10 optimization with different 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.005 0.020 0.023 0.026 0.029 0.032 0.035 0.038 0.040 0.050 ).
Due to the size of data, we only report suboptimal portfolios of 0-CDaR, 0.5-CVaR and standard deviation. We find that the optimal portfolio is consistently too conservative for all risk measures and input. The performance doesn't follow the expected trend that it would first monotonically increase until reaching the peak then monotonically decrease.
We report the performance in Table V and Fig.5 .
The portfolios have increasing return as well as risk from L1 to H4. The realized paths follow similar trend with rare cross. The optimal portfolios, as expected, are in the medium part of all paths. The 10-day 0.5-CVaR suboptimal portfolios have almost identical performance. Some paths are not visible due to overlap. For example, the constraints on return of H4 are sometimes unfeasible, leading to same performance with H3 portfolio.
We propose a Markov regime switching GARCH model with multivariate normal tempered stable innovation to generate simulated sample paths as input for portfolio optimization with CDaR, CVaR and standard deviation risk measure, aiming at modeling stylized facts, especially fat tail, accommodating regime switch and addressing tail risk in portfolio performance. We conduct in-sample and out-of-sample tests to examine the effectiveness of the model. In insample test, the model fits the residuals with high accuracy. In out-of-sample tests, we show Sharpe-like ratios, and that tail risk measures are slightly better than standard deviation. By combining MRS-MNTS-GARCH model and optimization with tail risk measures, the optimal portfolios successfully alleviate extreme left tail events. In a further study, we find that optimal portfolios tend to be conservative. Real optimal portfolios are those with higher constraint on return.