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

.


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.

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.

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.

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.

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 1 X and 2 X , 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.
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.
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, µ is the vector of means and ε is the vector of associated error.
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 i T is the target for the i th factor, Z is the vector of standardized variables, and i L is the vector of loadings composed of the correlations between the original variables and the ith factor.
Finally, the FMSE values for the i th factor can be calculated as shown in Equation 7, where i F is the value of the rotated factor score and i λ is the variance associated to each factor.
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.

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. Where i x indicates the proportion of the asset i in the portfolio, µ is the expected return of the asset, ij σ means the covariance between the returns of the i TH and j TH assets; 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 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).

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 , ϕ , and θ are the parameters of the model to estimate the mean t y , and the term ε is the error. On the other hand, in Equation 10, 0 ζ , i ζ , and λ are the parameters of the model to estimate the variance 2 σ , and ε is the error term. ( ) For: 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 corresponds to the score evaluated in t θ ; and t t θ Γ is defined, as shown in Equation 13.
The additional parameter γ is fixed and usually assumes values among the set (0, ½, 1). When

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.

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 1 2 , , , p w w w … are proportions of p components of a mixture design, then, Equation 14 can be written.
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. 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.

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: Where Y = n-dimensional response vector; X = n x p model matrix; and ε = aleatory error term. The variance of the parameters is given by Equation 17: And the covariance matrix is given by Equation 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.

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.

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 t P be the asset price and  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. We have also compiled some descriptive statistics about the series in order to better understand them. These statistics can be viewed in Table 1. 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.
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.

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

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.

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.
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.

( )
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 1 w , 2 w , 3 w , 4 w , and 5 w , respectively. 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.   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.
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, 1 F T and 2 F T , 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 1 λ e 2 λ are respectively equal to 1.97 and 1.00, as shown in Table 4. . 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.   Finally, the overlaid contour plots of FMSE 1 and FMSE 2 can be viewed in Figure 6, where the white area represents feasible solutions.

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.
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. 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.

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.