A multiobjective portfolio optimization for energy assets using D-Optimal design and mixture design of experiments

Abstract

Paper aims

Frequently, researchers try to find a better way to allocate assets in order to have maximum return and low variability in a portfolio as diverse as possible. This paper aims to apply D-Optimal Design in the context of Mixture Design and portfolio optimization to efficiently select the runs of the proposed experimental design.

Originality

A new approach to find the optimal weights that maximize the returns and minimize the risk using D-Optimal Design was used. A multi-response optimization problem considering returns, variability and entropy as functions of the weights was proposed. However, as there is a significant correlation between the objective functions, a Factor Analysis combined with FMSE to dimensionality reduction was used.

Research method

All the steps for both stages of the methodology applied in this paper are presented below: select real time series; predict one step ahead; generate a D-Optimal mixture design; apply weights and generate mathematical models; solve the optimization problem.

Main findings

Using the desirability method, the optimal values were determined, obtaining approximately 79% for the compound desirability function. The proposed method presented a 16.80% higher return with a 4.98% higher risk exposure if compared against Naïve method.

Implications for theory and practice

The proposed methodology can be applied to any portfolio optimization study. Mixture Design studies have already been proposed for modeling portfolio optimization problems. However, the D-Optimal Design proved to be adequate, which represents less computational effort.

Keywords
Portfolio optimization; D-Optimal design; Mixture design; Multivariate analysis

1. Introduction

Portfolio optimization is a term that can be understood as a problem whose main object is to encounter the most appropriate way to allocate the assets so that maximum returns and the smallest variability can be obtained (Garg & Deep, 2019). Markowitz first introduced the mean variance model in which the risk of the model is defined as the variance of the returns and nowadays many optimization portfolio studies are encountered in the literature (Çela et al., 2021; Li et al., 2019; Nguyen & Lo, 2012).

A commonly used criterion to choose the best portfolio is based on maximum expected return and minimum variance (Li & Shu, 2011; Mendes et al., 2016). However, many other approaches can be found. Entropic value-at-risk (EVaR) was applied in (Ahmadi-Javid & Fallah-Tafti, 2019) as a risk measure in sample-based portfolio optimization designing an efficient algorithm based on differentiable convex optimization.

According to (Mercurio et al., 2020), there are mainly five problems for the simplistic use of the method made by Markowitz: 1) optimal solutions attributing great weight to high-risk assets; 2) does not consider the structural dependence of the assets under analysis; 3) drastic variations in optimal solutions when adjusting the inputs of the model; 4) strong dependence that the returns present multivariate normal distribution; 5) difficulty in estimating the covariance matrix and expected returns.

A research developed by (Garg & Deep, 2019) demonstrated the effectiveness of using some variants of Biogeography Based Optimization to solve the mean variance model of portfolio optimization. The problem has great complexity since it presents a number of local optima. The data used in the study was extracted from Indian National Stock Exchange, Mumbai considering the period starting in April 2015 to the end of March 2016.

According to (Milhomem & Dantas, 2020), although optimization methods have been developed over the years, investors and academics often prefer deterministic techniques over the stochastic ones due to they are easier to understand and to apply.

A great opportunity in this area is to apply techniques commonly used to optimize manufacturing processes, such as mixture designs and the multiobjective optimization method called desirability as in (Rocha et al., 2017; Tillmann et al., 2010) to optimize portfolios composed of different assets. The study developed in (Mendes et al., 2016) aimed to optimize a portfolio consisting of distinct heteroskedastic assets of world crude oil spot prices using ARMA-GARCH models. After modeling the time series, the authors predicted the returns and the risk and used a mixture design to estimate different values of return and risk according to the weights specified in this design. Then, optimal weights values which optimizes the return and minimize the risk of the portfolio were obtained.

A portfolio optimization study, considering seven contract candidates from the Brazilian energy market, was performed in (de Oliveira et al., 2011). A novel approach using mixture design experiments to build up nonlinear models for the risk is presented. The desirability was used as the multi-response optimization technique in order to maximize the return and minimize the risk of the portfolio. An entropy objective function was added to ensure portfolio diversification, serving as an adaptation to the Markowitz model.

Advancing on the aforementioned works, this article proposes to apply an experimental design called D-Optimal Design. According to (Montgomery, 2017), the D-Optimal design applies a criterion for selecting the points of the experimental design in such a way that it results in the minimization of the confidence interval set for the regression coefficients.

According to (Rather et al., 2017) the construction of an optimal portfolio has to take into account a careful prediction of the variables considered for the construction of the model. This work therefore seeks to compare the prediction given between ARMA-GARCH models and GAS models, highlighting the importance of having accurate methods during portfolio optimization.

Generalized autoregressive conditional heteroskedasticity (GARCH) models, which is a natural generalization of autoregressive conditional heteroskedasticity (ARCH) (Bollerslev, 1986), are the most commonly used models and were also applied in (Chan & Grant, 2016) aiming to compare this class of models with stochastic volatility models using real series of oil, petroleum product and natural gas prices. Regarding the considered series, the authors concluded that stochastic volatility models outperform the traditional GARCH models in most cases.

On the other hand, another method used to model and predict a heteroskedastic time series is the generalized autoregressive score (GAS) model, a method whose main distinctive characteristic is that it uses a score function to drive time-variation with regards to non-linear models’ parameters (Ardia et al., 2019). GAS is applied in (Creal et al., 2013) as an applicable model to capture time variation in parameters. The authors have their studies based on dynamic copula models and multivariate marked point process models.

As could be observed, no previous study considered an analysis using D-Optimal Design in the context of portfolio optimization combined with mixture design and desirability, aiming to maximize the returns, minimize the risk and obtain a portfolio as diversified as possible. Moreover, this work compares GARCH and GAS in the selected series, which have heavy tails, through Rooling Window to define which method will be used to predict the risk and return of the selected series.

The remainder of this paper is organized as follows. Section 2 presents the materials and methodology used both to develop the analysis to optimize the considered portfolio using D-Optimal Design in the context of Mixture Design. In section 3, the literature review encompassing the contents which are necessary to understand the main concepts discussed in the article is presented. Section 4 presents the results and the discussions. Finally, section 5 presents the main conclusions of this work.

2. Materials and methods

The main purpose of the paper is to optimize a portfolio consisting of 5 distinct time series of assets in the energy sector, using D-Optimal design to support the proposed mixture design. Thus, it will be necessary to make predictions of these series, and it is important to apply an accurate method to measure both the mean and the variance. Afterwards, the best method will be applied during the optimization of the real portfolio.

Thus, the methodology applied in this scientific article consists of selecting real time series; to compare two suitable adjustment and forecasting methods for series with heavy tails; predict both return and risk one step ahead; generate a mixture design, through D-Optimal Design, in order to model the multi-objective optimization problem (maximizing return, minimizing risk and maximizing entropy) and solve the optimization problem through the technique of Desirability.

The classic Markowitz model of mean-variance, with the addition of the entropy function, was used for the modeling, since the present work seeks to apply D-Optimal Design to model the objective functions.

All the steps for both stages of the methodology applied in this paper are presented below. The next subsections present a detailed explanation of each step.

2.1. Select real time series

In this first step, n real time series, which define the portfolio to be optimized, must be selected. As risk and return forecasts are important for portfolio optimization, two models suitable to predict the variance (as a proxy for volatility) of a given time series were compared using a rolling window of size one.

2.2. Predict one step ahead

For each considered time series, the researcher predicts one step ahead for the returns values. The best method selected is used here to make the predictions.

2.3. Generate a mixture design

A mixture design is generated considering n factors. A simplex lattice design is recommended, since this design was successfully used in a time series optimization forecast problem in (Bacci et al., 2019). Nevertheless, it is convenient to define lower and upper bounds to the values of the mixture design as will be explained later.

2.4. Apply weights and generate mathematical models

Each run of the mixture design represents a set of m values of weights which multiplies the predictions performed in k. The expected return and variability of the portfolio will be calculated over a linear combination. Let X be the linear combination of two distinct variables $X_1$ and $X_2$, as shown in Equation 1. Then, the expected value and the variability of X can be calculated as shown in Equations 2 and 3, respectively. Hence it is possible to create a model for the expected value of the return and the variability as functions of the weights.

$X=a{X}_1+b{X}_2$ (1)

$E\left(X\right)=aE\left(X_1\right)+bE\left(X_2\right)$ (2)

$Var\left(X\right)=a^{2}Var\left(X_1\right)+b^{2}Var\left(X_2\right)+2abCov\left(X_1,X_2\right)$ (3)

Besides the values of return and variability, which are supposed to be maximized and minimized, respectively, we also have the entropy. According to (Mendes et al., 2016), the entropy measures the diversification of the portfolio, thus the greater the better and it is calculated as in Equation 4.

$S=-{\displaystyle \sum }_{i=1}^n{w}_{i}log\left(w_i\right)$ (4)

We highlight that this represents a multiobjective optimization problem and it is up to the researcher to decide which method is suitable for the purposes of the study. Nevertheless, a common issue in this context is to have correlated responses. In these situations, it is possible to apply multivariate techniques, such as factor analysis, in order to work with uncorrelated rotated factor scores. According to (Johnson & Wichern, 2007), the orthogonal factor model can be written as shown in Equation 5, where$X$ is the vector of the original variables, $F$ represent the vector rotated factor scores, $L$ is the matrix of loadings, $\mu$ is the vector of means and $\epsilon$ is the vector of associated error.

$X=\mu +LF+\epsilon$ (5)

The rotated factor scores can well represent the original correlated variables, however, in some cases, a single factor explains variables that present distinct optimization directions. For instance, return, to be maximized and variability, to be minimized. In order to overcome this problem, the technique of Factor Mean Square Error (FMSE) proposed by (Leite, 2019) as a variation of Multivariate Mean Square Error (MMSE) developed in (Paiva et al., 2009) can be applied as in (Luz et al., 2021). Initially, the targets for the factors are calculated as shown in Equation 6, where $T_i$ is the target for the i^th factor, Z is the vector of standardized variables, and $L_i$ is the vector of loadings composed of the correlations between the original variables and the ith factor.

$T_i=Z^T{L}_i$ (6)

Finally, the $FMSE$ values for the i^th factor can be calculated as shown in Equation 7, where $F_i$ is the value of the rotated factor score and ${\lambda }_i$ is the variance associated to each factor.

$FMS{E}_i={\left[F_i-T_i\right]}^2+{\lambda }_i$ (7)

2.5. Solve the optimization problem

In order to perform the second stage, we selected 5 real time series related to the energy sector, they are: West Texas Intermediate (WTI) in dollars per barrel, Europe Brent (Brent) in dollars per barrel, New York Harbor No. 2 Heating Oil (Heating) in dollars per gallon, Mont Belvieu, TX Propane (Propane) in dollars per gallon, and New York Harbor Conventional Gasoline Regular (NYHCGR) in dollars per gallon. They can all be found in U.S. Energy Information Administration (EIA) website. We highlight that all the series consider spot price FOB. Figure 1 shows the plots associated to each one of them.

Figure 1
Time series plots of the real time series of the energy sector.

3. Background and literature review

3.1. Portfolio optimization

Allocating resources in optimal portfolios is part of the study area of several scientific subjects. Markowitz's mean-variance model (Markowitz, 1952) can be interpreted in two concomitant ways. The investor wishes to maximize the return given a level of risk or given a desired return the investor wishes to minimize the risk of his stock portfolio. Markowitz's model considers absence of transaction costs. (Ross et al., 2017) define a portfolio being a set of assets, securities, real estate, investment projects, etc., whose objective is to reduce risk through a correct diversification of capital applied in different alternatives of financial capital allocation.

According to Anagnostopoulos & Mamanis (2010), considering a multi-objective approach, the problem can be described mathematically as Equation 8.

$\max \mu \left(x\right)={\displaystyle \sum }_{i=1}^n{w}_i{\mu }_i$

$\text{min}\rho \left(\text{x}\right)={\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^n{w}_i{w}_j{\sigma }_{ij}$ (8)

$St:$

${\displaystyle \sum }_{i=1}^n{w}_i=1$

${\displaystyle \sum }_{i=1}^n{\delta }_i\le K$

$l_i{\delta }_i\le x_i\le u_i{\delta }_i,i=\mathrm{1,}\mathrm{2,}\dots ,n$

${\delta }_i\in \left\{\mathrm{0,1}\right\},i=\mathrm{1,}\mathrm{2,}\dots ,n$

Where $x_i$ indicates the proportion of the asset i in the portfolio, $\mu$ is the expected return of the asset, ${\sigma }_{ij}$ means the covariance between the returns of the i^th and j^th assets; ${\delta }_i$ is a dummy variable that receives 1 if the asset is chosen for the portfolio and 0 if it is not, n is the number of assets included in the portfolio; K is the number of assets (or series) available. In the present work, the variable ${\delta }_i$ was not considered since all series were considered for the composition of the portfolio. In this way, the Markowitz model served as a basis for modeling the response surface given by the design of mixtures.

(Merton, 1971) works the optimization problem in relation to continuous time under uncertainty. One of the main advantages of the author's approach is that the problem in question can be analyzed through two stochastic processes: Brownian motion functions and Poisson processes, reducing the number of parameters to be estimated.

A possible solution to a problem of resource allocation in a portfolio of assets is the NAIVE solution, which corresponds to dividing the weight equally for n assets (1/n). According Behr et al. (2013) this strategy can be suitable in a context of great uncertainty, being a good strategy for investors who have a great aversion to risk.

According (Mansini et al., 2014), there is a family of medium-risk measures that can be used. One of the widely used risk measures was proposed by Uryasev (2000) which is called Conditional Value-At-Risk (CVaR). The CVaR corresponds to the weighted average of the extreme losses that can occur in a distribution of probability of returns.

Regarding the optimization methods, linear optimization methods with constraints can be used, such as the Simplex algorithm and interior point methods in order to reduce the computational effort. Other forms of optimization can be applied to optimization problems, such as heuristic methods (Cura, 2009; Hu et al., 2015).

3.2. Generalized autoregressive conditional heteroskedasticity models

In the financial literature the uncertainty is a crucial concept since the risk that quantifies them is presented in many models such as the Capital Asset Pricing Model (CAPM) and in the modern portfolio optimization theory and in general the notion of risk of an asset is associated with the historical volatility of its returns.

However, the volatility is a latent variable which is not observable and therefore must be computed. According to (Tsay, 2005) the unobservability of volatility makes it difficult to evaluate the forecasting performance of conditional heteroscedastic models and although not directly observable, presents some characteristics (stylized facts) that are commonly observed in asset returns, such as clustering and leveraging.

One model extensively found in the financial econometric literature and largely employed by practitioners is the Generalized autoregressive conditional heteroskedasticity (GARCH) model proposed by (Bollerslev, 1986).

The GARCH family models are observational models (Creal et al., 2013) and demonstrate high popularity due to the simplification of the evaluation of likelihood functions and possible prediction of parameters. In models directed to parameters, the parameters themselves are stochastic, with their own sources of error.

However, in many cases it is reasonable to assume that the variance depends on the variability of many previous periods, leading to the AR-GARCH model, which mathematically consists of Equations 7 and 8 of a GARCH model(l, s). In Equation 9, $c$, $\phi$, and $\theta$ are the parameters of the model to estimate the mean $y_t$, and the term $\epsilon$ is the error. On the other hand, in Equation 10, ${\zeta }_0$, ${\zeta }_i$, and $\lambda$ are the parameters of the model to estimate the variance ${\sigma }^2$, and $\epsilon$ is the error term.

$y_t=c+{\displaystyle \sum }_{i=1}^r{\phi }_i{y}_{t-i}+{\displaystyle \sum }_{j=1}^m{\theta }_i{\epsilon }_{t-j}+{\epsilon }_t$ (9)

For: ${\epsilon }_t={\sigma }_i{Z}_t;Z_t~N\left(\mathrm{0,}1\right)$.

$Var\left({\epsilon }_t\text{|}{\epsilon }_{t-1},\dots ,{\epsilon }_{t-j}\right)={\zeta }_0+{\displaystyle \sum }_{i=1}^p{\zeta }_i{\epsilon }_{t-j}^2+{\displaystyle \sum }_{i=1}^q{\lambda }_i{\sigma }_{t-j}^2$ (10)

For: ${\epsilon }_t~N\left(\mathrm{0,}{\sigma }_t^2\right)$.

3.3. Generalized autoregressive scores models

(Creal et al., 2013) and (Harvey & Sucarrat, 2014) argue that the models proposed to model series with conditional variances are difficult to estimate, besides not taking into account the form of conditional distribution of data.

(Creal et al., 2013) proposed, therefore, the use of conditional density function score as a guide for changing the variation-time of the parameters of the analyzed time series process. The authors indicate that a great advantage of this method is that the estimation by maximum likelihood is direct. An additional utility of the GAS model is its versatility, allowing a good adjustment of the modeling in series with several characteristics (Ardia et al., 2019).

Be $y_t\in {\Re }^N$ for a N-dimensional vector at time t, with a conditional distribution $y_t|y_{1:t-1}~p\left(y_t{\theta }_t\right)$. For $y_t|y_{1:t-1}\equiv {\left(y_1^T,\dots ,y_{T-1}^T\right)}^T$, and ${\theta }_t\in \Theta \subseteq {\Re }^j$ is a vector of parameters that vary in time and depends only of $y_{1:t-1}$ and another set of additional parameters ξ for the entire t. The ${\theta }_t$ parameters are updated by the score defined in Equation 11.

${\theta }_{t+1}=k+A{s}_t+B{\theta }_t$ (11)

Where$\kappa$, A, and B are matrices of coefficients and ${\theta }_t$ is a proportional vector to the score. Thus, $s_t$ can be defined as Equation 12.

$s_t\equiv S_t\left({\theta }_t\right){\nabla }_t\left(y_t,{\theta }_t\right)$ (12)

Where $S_t({\theta }_t)\equiv {\Gamma }_t{\theta }_t^{-\gamma }$ is a defined positive matrix of order J x J known at time t; ${\nabla }_t\left(y_t,{\theta }_t\right)\equiv \partial logp\left(y_t,{\theta }_t\right)/\partial {\theta }_t$ corresponds to the score evaluated in ${\theta }_t$; and ${\Gamma }_t{\theta }_t$ is defined, as shown in Equation 13.

${\Gamma }_t\left({\theta }_t\right)\equiv E_{t-1}\left[{\nabla }_t\left(y_t,{\theta }_t\right){\nabla }_t{\left(y_t,{\theta }_t\right)}^T\right]$ (13)

The additional parameter $\gamma$ is fixed and usually assumes values among the set (0, ½, 1). When $\gamma =0$, there is no scale and $S_t=I$, considering I an identity matrix of appropriate size. If $\gamma$ = 0 or $\gamma$ = ½ the conditional score ${\nabla }_t\left(y_t,{\theta }_t\right)$ is multiplied by the inverse of the covariance matrix ${\Gamma }_t\left({\theta }_t\right)$.

Reparameterization occurs when observing a new realization of the variable under analysis. In Equation 13 the ${\theta }_t$ is linear and, therefore, is unconstrained. However, considering, for example, that the scale vector follows a Student t-distribution, a non-negativity constraint is required. The parameters $\kappa$, A, and B are easily numerically estimated by maximum likelihood.

3.4. Design of experiments

Design of experiments (DOE) is an invaluable technique where the experimenter performs a certain number of designed experiments and it allows to make reliable and important statistical conclusions. The number of experiments to be performed is certainly smaller than those in a trial-and-error approach.

The 2-levels or $2^k$ factorial design, where k represents the number of factors being considered, is widely used in scientific applications. In this type of design, a first order model with interactions is generated, helping to understand the main effects of the factors over the objective function. Furthermore, it is also possible to understand how the factors jointly affect the response (Montgomery, 2017).

In a $2^k$ factorial design, the levels of the factors are usually defined as low and high. All the possible combinations of these levels are investigated in a series of tests (Yondo et al., 2018). In these tests, the levels of the considered factors are changed by the researcher in a controlled way in order to observe the effects of each factor and interactions (Staiculescu et al., 2005).

Usually, in a factorial design the points, representing each run of the structured design, are set at the vertices of a hypercube (Montgomery, 2017). It is worth mentioning that (1) indicated the point where all factors are set at their lowest levels and in all the other vertices the letters indicate which factors are set at their highest levels. According to (Yondo et al., 2018), industrial experiments have been focusing mainly on 2-levels and 3-levels designs, resulting in $2^k$ and $3^k$ experiments, respectively.

3.5. Mixture design

According to (Montgomery, 2017), the design of mixtures corresponds to the analysis of experiments in which the factors correspond to components or ingredients of a mixture, in which the levels are not independent. That is, if $w_1,w_2,\dots ,w_p$ are proportions of p components of a mixture design, then, Equation 14 can be written.

${\displaystyle \sum }_{i=1}^n{w}_i=1$ (14)

Comparing Equation 14 with the constraints of Equation 8, it is noted that portfolio optimization can be clearly treated as a mixture design problem. In an experiment analysis problem, explanatory variables (or factors) are fixed at certain levels to study their behavior in one or more response variables. In a mixture design, the proportions are modified to analyze the behavior of the response variable.

In a mixture problem with three components the region of feasible solutions and the coordinate system can be represented by a triangular plane. Therefore, we can estimate via Ordinary Least Squares a response surface designed to parallel cuts to the triangular base (Lawson, 2014). In general, it is arranged that there is a functional relationship between the response variables and the proportions, so that, Equation 15 can be written.

$y_i=f\left(w_1,w_2,\dots ,w_n\right)$ (15)

Where $y_i$ is the dependent variable; $w_i$, with $i=\mathrm{1,}\mathrm{2,}\dots ,n$, are the explanatory variables or proportions; and finally, $f$ is the functional form.

The mixing problem then consists of determining a model that represents the response surface. Considering the constraints imposed on the model, one can choose from the following standard experiments: simplex lattice, simplex centroid design, extreme vertices design.

3.6. D-optimal design

In recent years, optimal designs of experiments are being used to the detriment of classical techniques. There are many criteria for selecting optimal experimental designs. One can list: A-optimal, D-optimal, I-optimal and G-optimal criteria. The first two are criteria oriented to the estimation of parameters, while the last two are oriented to prediction.

While the A-Optimal design focuses on the variance of the estimated parameters, using the standard deviation to construct confidence intervals and apply hypothesis tests, the D-Optimal design minimizes the volume of multivariate confidence intervals (ellipses) around parameters (Montgomery, 2017).

Considering a general model exemplified by Equation 16, given by matrix notation:

$Y=X\beta +\epsilon$ (16)

Where$Y$= n-dimensional response vector; $X$ = n x p model matrix; and $\epsilon$= aleatory error term. The variance of the parameters is given by Equation 17:

$\widehat{\beta }={\left(X^{T}X\right)}^{-1}{X}^{T}Y$ (17)

And the covariance matrix is given by Equation 18:

${\sigma }^2{\left(X^{T}X\right)}^{-1}$ (18)

The inverse of the covariance matrix is Fisher's information matrix.

The runs are selected so that, the determinant of X'X is maximum. For first-order models, the full factorial corresponds to a D-Optimal design and the non-linearity presented in the modeling justifies the use of D-Optimal design. As this determinant is inversely proportional to the ellipsoidal confidence intervals, the D-optimal design minimizes the number of runs required compared to a full factorial design.

4. Results

This section presents the results obtained in each stage of the methodology previously explained. The decisions and results from each step of these stages are presented separately in order to make them easier to understand.

4.1. Select real time series

Since our goal is to predict the returns and volatilities of the considered times series, we choose to transform the series into their natural logarithms, which is more convenient to work with and provides some desirable properties as symmetry, monotonicity, scale invariance and time-additivity. The logarithmic transformation is in fact appropriate whenever the time series exhibits a standard deviation that increases linearly with the mean.

Therefore, let $P_t$ be the asset price and $R_t=\left(P_t-P_{t-1}\right)/P_{t-1}$ its ordinary net return. The log return (taking the natural logarithm) is computed, as shown in Equation 19:

$r_t=\ln \left(1+{\text{R}}_{\text{t}}\right)=\ln \left(\frac{P_t}{P_{t-1}}\right)=\text{ln}(p_t)-\text{ln}\left(p_{t-1}\right)$ (19)

According to Tsay (2005), the log-returns present some advantages over the ordinary ones since its statistical properties are more tractable. Figure 2 presents the series analyzed in log-return.

Figure 2
Time series plot of the returns of the considered time series.

We have also compiled some descriptive statistics about the series in order to better understand them. These statistics can be viewed in Table 1.

Table 1
Descriptive statistics of the energy sector time series.

Statistic	WTI	Brent	Heating	Propane	NYHCGR
Mean	-0.0866	-0.0417	-0.0397	-0.0079	-0.0356
Min.	-301.9661	-64.3699	-18.4602	-16.959	-29.9857
Max.	124.0941	41.2023	14.8625	19.9796	22.2221
s.d.	7.9679	3.1839	2.3554	2.9872	2.8679
Var.	63.4877	10.1370	5.5480	8.9231	8.2251
Kurtosis	1061.9193	111.7142	12.0888	9.0119	26.8329

Aiming to validate the result obtained in the first stage of the methodology, we modeled the real time series using both GAS and GARCH methods. Table 2 shows the results in terms of RMSE for the returns of these series.

Table 2
RMSE for the models of returns applying GARCH and GAS in real time series.

Time series	RMSE
	GARCH	GAS
1-WTI	17.1695	17.2824
2-Brent	5.9011	5.8902
3-Heating	3.2974	3.2857
4-Propane	4.0020	3.9441
5-NYHCGR	4.4551	3.4364

As can be observed in Table 2 the results are very close and both methods present a similar performance. However, GAS achieved better predictive performance in 4 of the 5 series analyzed. In view of this, GAS provides better results considering these series than GARCH, hence GAS was chosen to be applied in the next stage of this work.

4.2. Predict one step ahead

Using GAS we obtained Table 3 which contains the predictions for both returns and variability values.

Table 3
Predicted values of return and variance for the real time series.

Serie	Predicted return	Predicted variability
1-WTI	0.0104	2.8472
2-Brent	-0.0078	2.4902
3-Heating	-0.0131	2.2824
4-Propane	0.0527	3.7126
5-NYHCGR	0.0253	2.9456

4.3. Generate a mixture design

A simplex lattice mixture design with degree of lattice 10 was generated. Five components were considered resulting in 1006 runs. We established 0.05 and 0.80 as the lower and the upper bounds of the components, which guarantees the diversification of the portfolio, avoiding cases where a certain asset is not considered.

Nevertheless, in order to select the best points to generate the model of the returns, risk and entropy versus the values of weight, it was applied the D-optimal design technique. From the 1006 runs an optimal design was generated consisted of only 200 runs. Henceforth, it is the design that will be used in the next steps.

4.4. Apply weights and generate mathematical models

Each run of the mixture design was written as a linear combination of the predictions, where the coefficients or weights are defined by the mixture design. Then, it was possible to extract the expected value, and the variability through Equations 20 and 21, respectively.

$E_X={\displaystyle \sum \begin{array}{l}n\\ i\end{array}}={}_{1}w{}_i{x}_i$ (20)

where $w_i$ indicates the proportion to be invested in each asset represented by the time series.

$Var\left(E_X\right)={\displaystyle \sum }_{i=1}^n{w}_i{}^2{\sigma }_{X_i}^2+2{\displaystyle \sum }_{i,j:i<j}{w}_i{w}_{j}cov\left(X_i,X_j\right)$ (21)

We modeled the return, variability and the entropy as functions of the weights and these models can be viewed in Equations 22, 23 and 24.

$E\left(w\right)=0.01044w_1-0.007797w_2-0.01311w_3+0.05271w_4+0.02532w_5$ (22)

$Var\left(w\right)=2.82795w_1+2.46925w_2+2.22715w_3+3.52025w_4+2.90021w_5$ (23)

$$\begin{array}{c}S\left(w\right)=0.06711w_1+0.06640w_2+0.06673w_3+0.06680w_4+0.06642w_5\\ +1.4746w_1{w}_2+1.4746w_1{w}_3+1.4716w_1{w}_4+1.4737w_1{w}_5\\ +1.4772w_2{w}_3+1.4738w_2{w}_4+1.4779w_2{w}_5+1.4760w_3{w}_4\\ +1.4757w_3{w}_5+1.4749w_4{w}_5\end{array}$$ (24)

In order to better visualize the behavior of the models previously generated, the mixture contour plots and the mixture surface plots are presented in Figure 3, where A, B, C, D, and E represent the weights $w_1$,$w_2$, $w_3$, $w_4$, and $w_5$, respectively.

Figure 3
Mixture contour plot and mixture surface plot of the return, variability and entropy.

After performing a correlation analysis, it was possible to state that the values of returns and variability are extremely correlated, as depicted in Figure 4. In view of this, a factor analysis process is adequate to reduce the dimensionality of the problem, but mainly, in order to work with uncorrelated factor scores.

Figure 4
Correlation matrix for the investigated response.

Table 4 presents the sorted rotated factor loadings and communalities for the factor analysis performed considering the returns, the variance and the entropy. In addition, it was possible to observe that 99.4% of the variance of the dataset can be explained with only 2 factors.

Table 4
Loadings and communalities from factor analysis.

Variable	Factor 1	Factor 2	Communalities
E(w)	0.995	-0.007	0.991
Var(w)	0.994	0.047	0.991
S(w)	-0.02	-1.000	1.000
Variance	1.980	1.002	2.981
% Variance	0.660	0.334	0.994

The uncorrelated factor scores, composed of 200 observations, were stored and they will be used in the next analysis, instead of the original variables. However, Factor 1 comprises both responses returns and variability, thus, it explains a response to be maximized and a response to be minimized. In order to overcome this issue, the FMSE technique, as previously explained, was applied. The target for both factors, $T_{F_1}$ and $T_{F_2}$, were calculated using Equation 7 and they are equal to 0.62 and -1.90, respectively. Afterwards, the FMSE values were calculate for each one of the 200 observations of the dataset of unrotated factor scores, where ${\lambda }_1$ e ${\lambda }_2$ are respectively equal to 1.97 and 1.00, as shown in Table 4.

Hence, the new responses, FMSE values were modeled according to the mixture design, as shown in Equation 25 and Equation 26. These equations will be considered in the multiple response optimization problem in the next step. Moreover, the contour and surface plots were also obtained for the FMSE models are shown in Figure 5.

$\begin{array}{c}FMS{E}_1=1.221w_1+5.715w_2+9.601w_3+8.176w_4+0.217w_5-4.913w1w2\\ -7.557w_1{w}_3-13.925w_1{w}_4-2.178w_1{w}_5-2.068w_2{w}_3\\ -26.714w_2{w}_4-6.865w_2{w}_5-34.262w_3{w}_4-10.949w_3{w}_5\\ -8.353w_4{w}_5\end{array}$ (25)

$R^2=98.94\%$

$R_{adj}^2=98.85\%$

$R_{pred}^2=98.58\%$

Figure 5
Mixture contour and surface plots for the FMSE₁ and FMSE_2.

$\begin{array}{c}FMS{E}_2=31.121w_1+31.047w_2+31.063w_3\\ +30.313w_4+30.419w_5\\ -87.05w1w2-86.94w_1{w}_3\\ -86.37w_1{w}_4-86.12w_1{w}_5\\ -85.92w_2{w}_3-86.04w_2{w}_4\\ -85.55w_2{w}_5-86.07w_3{w}_4\\ -85.84w_3{w}_5-85.03w_4{w}_5\end{array}$ (26)

$R^2=95.15\%$

$R_{adj}^2=94.78\%$

$R_{pred}^2=94.11\%$

Finally, the overlaid contour plots of FMSE₁ and FMSE₂ can be viewed in Figure 6, where the white area represents feasible solutions.

Figure 6
Overlaid contour plots of FMSE1 and FMSE2.

4.5. Optimization

The multiobjective optimization problem is now based on minimizing the FMSE values. In the desirability method it is necessary to establish some parameters which are low and high limits, targets, weights and importance to each objective function.

As far as the upper limits goes, they were set as the mean value of each response, whereas the targets and the lower limit were defined as zero. These settings were defined considering the response values calculate in the 200 runs of the MDE. Equal weights and importance were established to the three objective functions. Table 5 summarizes these settings. It is important to mention that the weight values indicate how the algorithm will emphasize or not the importance of hitting the target. Thus, a weight of 0.1 indicates that the algorithm will place less emphasis on the target. Furthermore, the importance values define the effect of each response over the composite desirability, which is the general function being optimized. Setting all the values of importance as 1, all the objective functions will have the same effect.

Table 5
Desirability settings.

Response	Goal	Lower	Target	Upper	Weight	Importance
FMSE1	Minimize	0	0	1.4	0.1	1
FMSE2	Minimize	0	0	4.6	0.1	1

After solving the optimization problem, it was obtained a composite desirability of 100% as can be viewed in Figure 7. The values in red represents the optimal weights to be applied on the time series assets.

Figure 7
Optimal values for the weights after using the desirability method.

Considering these weights, the values for return, variability and entropy will be 0.02, 2.92, and 0.66, respectively.

In order to create a comparison with other methods used, the proposed method was compared with the naive method of asset portfolio allocation. As exposed, as the proposed assets have high variability, the naive method is indicated for these cases.

As there are five assets under analysis, the Naive method proposes the allocation of 20% of the resources in each asset. The method proposed in this article suggests the allocation of 23.18% in the WTI asset, 19.01% in the Brent asset, 5.77% in the Heating asset and 26.02% in the Propane and NYHCGR assets. Table 6 summarizes this information.

Table 6
Allocation of resources via the NAIVE strategy and the one proposed by this paper.

Asset	NAIVE	D-Optimal Design
WTI	20%	23.18%
Brent	20%	19.01%
Heating	20%	5.77%
Propane	20%	26.02%
NYHCGR	20%	26.02%

Considering the time t+1 to t+30, the return generated by the Naive strategy is 1.29%, while that of the model proposed by this article was 1.51%. As for the variability, given by the estimate of the GAS method, the variability of the portfolio composed by the Naive method was 2.8511, while that of the portfolio proposed by this work was 2.9930. In other words, the method proposed by this article achieved a 16.80% higher profitability, for an increase in risk exposure of 4.98% in relation to the NAIVE strategy.

5. Conclusions

The presented paper aimed to perform an optimization portfolio study whose objective was maximize the return, minimize the risk (variability), in a portfolio as diversified as possible.

Due to the importance of making an accurate forecast, initially a comparative study on the series selected between GARCH and GAS methods was applied, applying Rolling Window technique, to determine which of the methods to use in predicting the return and risk of the selected series.

Since the GAS model performed better compared to the GARCH model, GAS method was applied to make one step ahead predictions of the real time series. The simplex lattice design with degree of lattice equal to 10, which is a type of mixture design, was considered to generate different sets of weights (each run of the design) used to weight the one step ahead predictions.

Nevertheless, a total of 1,006 runs are supposed to be considered in this case. Thus, the D-optimal design proved to be an adequate technique in this sense, because only 200 runs were considered, providing a satisfactory understanding about how the weight influenced the final results. Thus, it was possible to develop models to the return, variability, and entropy as functions of these weights.

As the objective functions showed high correlation, a Factor Analysis with FMSE to dimensionality reduction was applied. The desirability multi-response optimization method was applied on the reduced functions, in order to obtain the maximum return, minimum variability and maximum entropy. The composite desirability was equal to 1, which means that the objective was achieved and values equal to 0.02, 2.92 and 0.66 were obtained for the return, variability and entropy, respectively. As suggestions for new works, other optimal design models can be applied (such as A-Optimal and I-Optimal), in addition to the application of backtesting for the practical validation of the proposed theoretical model.

Finally, we compare the capital allocation in the analyzed assets between the method proposed by this paper with the equitable allocation that is suitable for investors with a moderate profile with data from one month after the analyzed period. The proposed method presented a 16.80% higher return with a 4.98% higher risk exposure.

Acknowledgements

The authors would like to thank the Brazilian agencies of CAPES, CNPq and FAPEMIG for supporting this research.

• How to cite this article: Leal, G. S., Romão, E. L., Reis, D. L. P. E., Balestrassi, P. P., & Paiva, A. P. (2022). A multiobjective portfolio optimization for energy assets using D-Optimal design and mixture design of experiments. Production, 32, e20210119. https://doi.org/10.1590/0103-6513.20210119

References

• Ahmadi-Javid, A., & Fallah-Tafti, M. (2019). Portfolio optimization with entropic value-at-risk. European Journal of Operational Research, 279(1), 225-241. http://dx.doi.org/10.1016/j.ejor.2019.02.007
  » http://dx.doi.org/10.1016/j.ejor.2019.02.007
• Anagnostopoulos, K. P., & Mamanis, G. (2010). A portfolio optimization model with three objectives and discrete variables. Computers & Operations Research, 37(7), 1285-1297. http://dx.doi.org/10.1016/j.cor.2009.09.009
  » http://dx.doi.org/10.1016/j.cor.2009.09.009
• Ardia, D., Boudt, K., & Catania, L. (2019). Generalized autoregressive score models in R: the GAS package. Journal of Statistical Software, 88(1). http://dx.doi.org/10.18637/jss.v088.i06
  » http://dx.doi.org/10.18637/jss.v088.i06
• Bacci, L. A., Mello, L. G., Incerti, T., Paiva, A. P., & Balestrassi, P. P. (2019). Optimization of combined time series methods to forecast the demand for coffee in Brazil: a new approach using Normal Boundary Intersection coupled with mixture designs of experiments and rotated factor scores. International Journal of Production Economics, 212(January), 186-211. http://dx.doi.org/10.1016/j.ijpe.2019.03.001
  » http://dx.doi.org/10.1016/j.ijpe.2019.03.001
• Behr, P., Guettler, A., & Miebs, F. (2013). On portfolio optimization: imposing the right constraints. Journal of Banking & Finance, 37(4), 1232-1242. http://dx.doi.org/10.1016/j.jbankfin.2012.11.020
  » http://dx.doi.org/10.1016/j.jbankfin.2012.11.020
• Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307-327. http://dx.doi.org/10.1016/0304-4076(86)90063-1
  » http://dx.doi.org/10.1016/0304-4076(86)90063-1
• Çela, E., Hafner, S., Mestel, R., & Pferschy, U. (2021). Mean-variance portfolio optimization based on ordinal information. Journal of Banking & Finance, 122, 105989. http://dx.doi.org/10.1016/j.jbankfin.2020.105989
  » http://dx.doi.org/10.1016/j.jbankfin.2020.105989
• Chan, J. C. C., & Grant, A. L. (2016). Modeling energy price dynamics: GARCH versus stochastic volatility. Energy Economics, 54, 182-189. http://dx.doi.org/10.1016/j.eneco.2015.12.003
  » http://dx.doi.org/10.1016/j.eneco.2015.12.003
• Creal, D., Koopman, S. J., & Lucas, A. (2013). Generalized autoregressive score models with applications. Journal of Applied Econometrics, 28(5), 777-795. http://dx.doi.org/10.1002/jae.1279
  » http://dx.doi.org/10.1002/jae.1279
• Cura, T. (2009). Particle swarm optimization approach to portfolio optimization. Nonlinear Analysis Real World Applications, 10(4), 2396-2406. http://dx.doi.org/10.1016/j.nonrwa.2008.04.023
  » http://dx.doi.org/10.1016/j.nonrwa.2008.04.023
• Garg, V., & Deep, K. (2019). Portfolio optimization using Laplacian biogeography based optimization. Opsearch, 56(4), 1117-1141. http://dx.doi.org/10.1007/s12597-019-00400-4
  » http://dx.doi.org/10.1007/s12597-019-00400-4
• Harvey, A., & Sucarrat, G. (2014). EGARCH models with fat tails, skewness and leverage. Computational Statistics & Data Analysis, 76, 320-338. http://dx.doi.org/10.1016/j.csda.2013.09.022
  » http://dx.doi.org/10.1016/j.csda.2013.09.022
• Hu, Y., Liu, K., Zhang, X., Su, L., Ngai, E. W. T., & Liu, M. (2015). Application of evolutionary computation for rule discovery in stock algorithmic trading: a literature review. Applied Soft Computing, 36, 534-551. http://dx.doi.org/10.1016/j.asoc.2015.07.008
  » http://dx.doi.org/10.1016/j.asoc.2015.07.008
• Johnson, R. A., & Wichern, D. W. (2007). Applied multivariate statistical analysis (6th ed.). New Jersey: Pearson.
• Lawson, J. (2014). Design and analysis of experiments with R New York: CRC Press. http://dx.doi.org/10.1201/b17883
  » http://dx.doi.org/10.1201/b17883
• Leite, R. R. (2019). Método de interseção normal à fronteira para modelos quadráticos de escores fatoriais rotacionais (Dissertação de mestrado). Universidade Federal de Itajubá, Itajubá.
• Li, Y., Jiang, X. F., Tian, Y., Li, S. P., & Zheng, B. (2019). Portfolio optimization based on network topology. Physica A, 515, 671-681. http://dx.doi.org/10.1016/j.physa.2018.10.014
  » http://dx.doi.org/10.1016/j.physa.2018.10.014
• Li, Z., & Shu, H. (2011). Optimal portfolio selection with liability management and Markov switching under constrained variance. Computers & Mathematics with Applications, 61(8), 2271-2277. http://dx.doi.org/10.1016/j.camwa.2010.09.045
  » http://dx.doi.org/10.1016/j.camwa.2010.09.045
• Luz, E. R., Romão, E. L., Streitenberger, S. C., Gomes, J. H. F., Paiva, A. P., & Balestrassi, P. P. (2021). A new multiobjective optimization with elliptical constraints approach for nonlinear models implemented in a stainless steel cladding process. International Journal of Advanced Manufacturing Technology, 113(5-6), 1469-1484. http://dx.doi.org/10.1007/s00170-020-06581-3
  » http://dx.doi.org/10.1007/s00170-020-06581-3
• Mansini, R., Ogryczak, W., & Speranza, M. G. (2014). Twenty years of linear programming based portfolio optimization. European Journal of Operational Research, 234(2), 518-535. http://dx.doi.org/10.1016/j.ejor.2013.08.035
  » http://dx.doi.org/10.1016/j.ejor.2013.08.035
• Markowitz, H. (1952). Portfolio selection. The Journal of Finance, 7(1), 77-91.
• Mendes, R. R. A., Paiva, A. P., Peruchi, R. S., Balestrassi, P. P., Leme, R. C., & Silva, M. B. (2016). Multiobjective portfolio optimization of ARMA-GARCH time series based on experimental designs. Computers & Operations Research, 66, 434-444. http://dx.doi.org/10.1016/j.cor.2015.05.001
  » http://dx.doi.org/10.1016/j.cor.2015.05.001
• Mercurio, P. J., Wu, Y., & Xie, H. (2020). An entropy-based approach to portfolio optimization. Entropy (Basel, Switzerland), 22(3), 1-17. http://dx.doi.org/10.3390/e22030332 PMid:33286106.
  » http://dx.doi.org/10.3390/e22030332
• Merton, R. C. (1971). Optimum consumption and portfolio rules in a continuous-time model. In W. T. Ziemba & R. G. Vickson (Eds.), Stochastic optimization models in finance (pp. 621–661.) New York: Academic Press.
• Milhomem, D. A., & Dantas, M. J. P. (2020). Analysis of new approaches used in portfolio optimization: a systematic literature review. Production, 30(1), 1-16. http://dx.doi.org/10.1590/0103-6513.20190144
  » http://dx.doi.org/10.1590/0103-6513.20190144
• Montgomery, D. C. (2017). Design and analysis of experiments (9th ed.). Hoboken: John Wiley & Sons.
• Nguyen, T. D., & Lo, A. W. (2012). Robust ranking and portfolio optimization. European Journal of Operational Research, 221(2), 407-416. http://dx.doi.org/10.1016/j.ejor.2012.03.023
  » http://dx.doi.org/10.1016/j.ejor.2012.03.023
• Oliveira, F. A., Paiva, A. P., Lima, J. W. M., Balestrassi, P. P., & Mendes, R. R. A. (2011). Portfolio optimization using Mixture Design of Experiments: Scheduling trades within electricity markets. Energy Economics, 33(1), 24-32. http://dx.doi.org/10.1016/j.eneco.2010.09.008
  » http://dx.doi.org/10.1016/j.eneco.2010.09.008
• Paiva, A. P., Paiva, E. J., Ferreira, J. R., Balestrassi, P. P., & Costa, S. C. (2009). A multivariate mean square error optimization of AISI 52100 hardened steel turning. International Journal of Advanced Manufacturing Technology, 43(7–8), 631-643. http://dx.doi.org/10.1007/s00170-008-1745-5
  » http://dx.doi.org/10.1007/s00170-008-1745-5
• Rather, A. M., Sastry, V. N., & Agarwal, A. (2017). Stock market prediction and Portfolio selection models: a survey. OPSEARCH, 54(3), 558-579. http://dx.doi.org/10.1007/s12597-016-0289-y
  » http://dx.doi.org/10.1007/s12597-016-0289-y
• Rocha, L. C. S., Paiva, A. P., Rotela Junior, P., Balestrassi, P. P., Campos, P. H. S., & Davim, J. P. (2017). Robust weighting applied to optimization of AISI H13 hardened-steel turning process with ceramic wiper tool: a diversity-based approach. Precision Engineering, 50, 235-247. http://dx.doi.org/10.1016/j.precisioneng.2017.05.011
  » http://dx.doi.org/10.1016/j.precisioneng.2017.05.011
• Ross, S. A., Westerfield, R., Jaffe, J., & Jordan, B. (2017). Essentials of corporate finance (9th ed). New York: McGraw-Hill Education.
• Staiculescu, D., Bushyager, N., Obatoyinbo, A., Martin, L. J., & Tentzeris, M. M. (2005). Design and optimization of 3-D compact stripline and microstrip Bluetooth/WLAN balun architectures using the design of experiments technique. IEEE Transactions on Antennas and Propagation, 53(5), 1805-1812. http://dx.doi.org/10.1109/TAP.2005.846820
  » http://dx.doi.org/10.1109/TAP.2005.846820
• Tillmann, W., Vogli, E., Baumann, I., Kopp, G., & Weihs, C. (2010). Desirability-based multi-criteria optimization of HVOF spray experiments to manufacture fine structured wear-resistant 75Cr 3C 2-25(NiCr20) coatings. Journal of Thermal Spray Technology, 19(1–2), 392-408. http://dx.doi.org/10.1007/s11666-009-9383-5
  » http://dx.doi.org/10.1007/s11666-009-9383-5
• Tsay, R. S. (2005).Analysis of financial time series New York: John Wiley & Sons. http://dx.doi.org/10.1002/0471746193
  » http://dx.doi.org/10.1002/0471746193
• Uryasev, S. (2000). Conditional value-at-risk: optimization algorithms and applications. In Proceedings of the IEEE/IAFE/INFORMS 2000 Conference on Computational Intelligence for Financial Engineering (CIFEr) (pp. 49-57). New York: IEEE.
• Yondo, R., Andrés, E., & Valero, E. (2018). A review on design of experiments and surrogate models in aircraft real-time and many-query aerodynamic analyses. Progress in Aerospace Sciences, 96, 23-61. http://dx.doi.org/10.1016/j.paerosci.2017.11.003
  » http://dx.doi.org/10.1016/j.paerosci.2017.11.003

Publication Dates

• Publication in this collection
  12 Sept 2022
• Date of issue
  2022

History

• Received
  01 Oct 2021
• Accepted
  05 Aug 2022