Forecasting using predictor selection from a large set of highly correlated variables

Abstract

The potential of correlation-based feature selection has been explored in selecting an optimal subset from a set of highly correlated predictors. This problem occurs, for example, in time series forecasting of economic indicators using regression models on multiple lags of a large number of candidate leading indicators. Greedy algorithms (forward selection and backward elimination) in such cases fail. To obtain the globally optimal solution, the feature selection problem is formulated as a mixed integer programming problem. To solve it, we use the binary cut-and-branch method. The results of simulation studies demonstrate the advantage of using the binary cut-and-branch method in comparison with heuristic search algorithms. The real example of the selection of leading indicators of consumer price index growth shows the acceptability of using the correlation-based feature selection method.

1. Introduction

Big data analytics includes the feature selection task [1, 2] for predictive modelling [3]. In many practical applications, candidate predictors correlate strongly. An example is the task of time series forecasting using leading indicators [4].

Lagged predictors are highly correlated. Fast and scalable univariate feature selection methods are not suitable under a given situation. They evaluate features individually, so the final subset includes many redundant strongly correlated features.

Multivariate methods take into account feature dependencies and try to discard not only irrelevant variables (which do not affect the response), but also redundant ones. Most often, the predictors are selected simultaneously with the construction of predictive models using embedded methods such as LASSO regression [4, 5]. It provides a sparse solution that includes only relevant features, which, however, is very sensitive to the regularization parameter.

In addition, stepwise regression is often used in time series forecasting [6]. It refers to the so-called “wrapper” methods. They select the optimal subset of features from all possible candidates simultaneously with the model estimation. Generally, this problem has exponential complexity in the number of features. In practice, to solve it, search approaches use greedy algorithms [7]. However, they do not usually produce an optimal solution, but approximate a globally optimal solution in a reasonable amount of time.

Finally, filter methods select variables regardless of the model. A Correlation-based Feature Selection is a well-known multivariate filter algorithm [8]. This approach is proposed for solving classification problems. Its applicability for predictor selection with a highly correlation of candidate predictors, in particular when selecting leading indicators, is poorly studied. This is the gap that we will attempt to address in our article. For this purpose, we first transform a correlation-based heuristic evaluation function optimization problem into a mixed integer programming problem. But in this problem the number of variables and constraints depends on the square of the number of features. With a branch and bound algorithm, the amount of computation becomes large. Therefore it is proposed to use the previously developed binary cut-and-branch method.

2. Correlation-based Feature Selection

Correlation-based Feature Selection (CFS) ranks feature subsets according to a correlation-based heuristic evaluation function [8]. The best subset contains predictors highly correlated with the response, yet uncorrelated to each other. Thus, the problem of feature selection is formulated as the following optimization problem:  Ri iSk  max ,

Sk k  2  rij

i, jSk ,i j where Ri is an absolute value of correlation coefficient between the response and the i-th feature, rij is an absolute value of correlation coefficient between the i-th and the j-th features, Sk is a subset of k features.

As for the time series of economic indicators, both the response and the predictors are usually quantitative. Therefore the Pearson product-moment correlation coefficient is applicable.

We reformulate the problem ( 1 ) as a problem of nonlinear integer programming: n  Ri2 xi  2 Ri Rj xi x j i1 i j n  xi  2 rij xi x j i1 i j where xi 0,1,i  1,..., n , n is the number of features. If xi  1 then the optimal subset contains the i max , x1, ,xn th variable, and xi  0 , otherwise.

The problem ( 2 ) is a polynomial fractional programming problem. Based on the transformation of feature selection problem proposed in [9], we replace the denominator in ( 2 ) by a positive continuous variable u. It leads to the equivalent polynomial problem. In addition, we convert a maximization problem into a minimization one. Thus, the problem is represented as follows:

n  Ri2 xiu  2 Ri Rj xi x ju  min i1 i j x1,...,xn ,u ( 1 ) ( 2 ) n  xiu  2 rij xi x ju  1 , i1 i j

u  0, xi 0,1.

Based on a linearization technique proposed in [10] to transform the terms xiu , xi x ju , we introduce variables zi , i  1,..., n , vij , i  1,..., n , j  1,..., n , i  j . Then we obtain the following mixed integer linear programming problem:

n  Ri2 zi  2 Ri Rjvij  ( 3 ) i1 i j min x1,...,xn ,u, z1,...,zn ,v11,...,vnn zi  0, vij  0, u  0, xi 0,1 , n  zi  2 rijvij  1, i1 i j

M  xi  1  u  zi  M 1 xi   u, zi  Mxi ,

M  xi  x j  2  u  vij  M  2  xi  x j   u, vij  Mxi ,vij  Mx j , where M is a large positive value.

Finally, the initial nonlinear problem ( 1 ) is reduced to a high-dimensional linear programming (LP) problem. The number of the new continuous variables is (n2+n+2)/2. They are added to the initial n binary variables. The number of constraints also depends on n2 and is (2n2+n+2).

3. Binary cut-and-branch method

The binary cut and branch method (BCBM) was originally developed to solve LP problems with Boolean variables [11] and then extended to the case of the General linear programming problem with mixed variables (milp) [12]. Any such problem, the special case of which is ( 3 ), can be represented as:  ( x)  c1T x  c2T y  const  max ,

A1x  A2 y  b , 0  x  1 , y  0 , ( 4 ) ( 5 ) ( 7 ) ( 8 ) ( 9 ) that the solution components x belong to one of the vertices of the unit hypercube of dimension n1 ; c1, x, 0, 1 are vectors of the same dimension; 0 is a zero vector; 1 is a vector of ones; const is a constant; and c1  0 . The vectors c2 , y, 0 have dimension n2 . Condition x  Ikn indicates that x belongs to the set of vectors of dimension n , the elements of which take integer values from the range [0  k 1] . Conformity of statements ( 4-6 ) and ( 3 ): Boolean variables x  I2n1 have the same meaning in both statements, continuous variables y  0 have the meaning of variables zi  0, vij  0, u  0 in statement ( 3 ). Matrices A1 and A2 in the constraints ( 4 ) are formed from the coefficients of the constraints in the problem ( 3 ).

In fact, any milp problem and a considerable part of mip problems can be compactly reduced to ( 1 )–( 3 ); see, e.g., [11].

Suppose x0, y0 is the solution of the relaxed problem ( 4 )–( 5 );  is the integer part of number; and  0   T x0 , where  j  0, 1 , j  1, n . Then any inequality of the form

 T x   0 ,  j  0, 1 , j  1, n ,  0   0  ,  0   T x0 , is called a binary cut (BC) for problem ( 4 )–( 5 ).

If x0 is part of the solution x0, y0 of the relaxed problem ( 4 )–( 5 ), then

 T x  0 , 0   T x0 .

Relation ( 8 ) can be a generating inequality if  j   iaij , i  0 , where aij ,i  I B are the iI B coefficients of the basis part A1 and i are the weights of the basis constraints. Specifically, if i are dual estimates for constraints ( 5 ), then  j  c j , j  1, n . A complementary system of BCs to constraints ( 5 ) is defined as

1D x   , where 1D   ij mD n1 , ij { 0, 1 }, j  1, n is the coefficient matrix of the complementary system and the vector  composed of the right-hand part of the constraints is defined from ( 7 ).

There are several ways to find out whether the BCs are valid [9,10]. Specifically, the binary cut and branch algorithm (BCBA) uses the following feature. We now define  T x   ( x0 )  1 , ( 10 ) where  ( x0 )   T x0  and x0 is part of the optimal solution x0, y0 of problem ( 4 ), ( 5 ), and ( 9 ).

 

If problem ( 4 ), ( 5 ), and ( 10 ) has a solution, then the cut  T x   ( x0 ) is invalid. Contrary wise, if problem ( 4 ), ( 5 ), and ( 10 ) has no solution due to the conflicting conditions ( 5 ) and ( 10 ), then  T x   ( x0 ) is a valid BC. This is the underlying feature of the BC synthesis procedure called selection in a set of the nearest cuts (SSNC) [11, 12]. We now describe this procedure.

We define an inequality ensuing from the basis system ( 5 ) and ( 9 ) through a permutation by arranging  in a descending order (denoted by  ). We consider a totality of n1 vectors, of dimension n1 : 1  1  (1,0,...,0) , ...,  j  (1,1,...,1,0,0,...,0) ( j original ones), ...,  n  (1,1,...,1) .

Each  j is set in correspondence with the value cs( j )  , j  1, n1 .  T j  2  j 2

The discrete function cs( j ) has a strict maximum and uniquely defines the priority of each of the alternative cuts with the coefficients  j . Adding the entire totality of these BCs to ( 9 ) and solving ( 4 ), ( 5 ), and ( 10 ), we can find out conflicting conditions (if there are any) to identify valid cuts. Then, if there are valid BCs, we select a single cut with the maximum value of cs( j ) . If there are no valid cuts, we select a BC that corresponds to the maximum of cs( j ) , j  1, n1 .

Another important feature of BCs is their radicality measure, which characterizes the depth of a cut of a given type. For a BC  T x   0 ,  j 0,1, j  1, n ,  0   0  ,  0   T x0 , we define the radicality r as the number of vertices of the unit hypercube cut off by the BC (BC system), assuming that the cut is valid.

In the general case, aT x  b ,

1 1 b  Ik1 , x  I2n , a  I2n , n1  a j x j  b, x j  I21, j  1, n1 , where a j 0,1} are the cut-off coefficients. j1

n1 For an arbitrary b  I1k , k   a j ,1  k  n1 , and Сkl  j1

k ! l!(k  l)! , we define

n1 k   a j ,1  k  n1

j1 rkb  2n1 k k  Сkl  max l b1

Relation ( 11 ) considers the unit hypercube vertices lying above the level aT x  b , i.e., belonging to the hyperplanes aT x  l with the right-hand parts ( b 1,b  2,...,k ). The aT x  b hyperplane itself contains 2n1 k Ckb vertices. The maximum radical BC is derived from ( 11 ) and  j  1, j  1, n with the possible exclusion of the minimum order relative to  if the sum of the coefficients in the lefthand part of the BC aT x  b without this exclusion is an integer number.

Regardless of which measure—closeness to the generating inequality or radicality is considered as a priority, the BCBA is as follows. or ( 11 ) 3.1. Binary Cut-and-Branch Algorithm

1. Suppose that we have obtained the solution of the original relaxed problem ( 4 ), ( 5 ), and ( 7 ): x0, y0 and  ( x0, y0 ) . If x0 are integers, the algorithm stops. Otherwise, it goes to step 2. stops. If the vertex with the maximum estimate  ( xq, yq ) contains integer xq , the solution ( xq , yq ) is the optimal one. The algorithm stops. Otherwise:

3. We create two new candidates for each of which we supplement the current matrix 1D for the step q with BCs ( 7 ) and ( 10 ) by the cut selection procedures (by the value of cs( j ) or by radicality ( 11 )): ( (t1) )T x   ( xq ) и ( (t1) )T x   ( xq ) 1, respectively.

4. We solve a pair of alternative subproblems with the cuts ( (t1) )T x   ( xq ) and ( (t1) )T x   ( xq ) 1.

5. We save their solution components xt 1 and x t 1 and the estimates  ( xt1, yt1) and  ( xt1, yt1) by adding them to a list of the tree vertices. If any of the candidates has no solution, it is withdrawn from the list of the vertices.

6. We increase the step number ( t : t 1 ) and go to step 2.

4. Simulation study results

The applicability of the CFS method to variable selection with a highly correlation of candidate predictors was investigated using the following model example.

The relevant features x1(m) , x2(m) , x3(m) , x5(m) , x6(m) were modeled as independent random variables with a standard normal distribution. The relevant features x4(m) , x7(m) were computed as follows: x(m)  x3(m)  e1, x7(m)  x(m)  x6(m)  e2 , 4 5 where e1, e2 are independently standard normally distributed.

The response was defined as

the redundant variables correlated with the relevant one, x1(r) ,..., x1(4r) ;  the irrelevant features 1,...,5 , 1,..., 4 .

A total number of variables was thirty. Simulation studies were repeated 1,000 times. The mixed integer linear programming problem contained 496 variables and 1,832 constraints. In order to solve the problem, greedy algorithms were used as an alternative: forward selection and backward elimination [13]. Their implementation in the R environment was used: the forward.search and backward.search functions provided by the FSelector package. The variable subset was evaluated by the objective function from ( 1 ).

Table 1 presents the results of simulation studies, the proportion of cases in which each variable is included in the subset. The irrelevant features were not included in the subset in any experiment. The following notation is used: forward — the forward selection, backward — the backward elimination, BCBM — the binary cut-and-branch method. forward

Method backward 0.927 0.929 0.962 0.98 0.977 0.969 0.999 0.057 0.066 0.053

Variable x(r) 4 x(r) 5 x(r) 6 x(r) 7 x(r) 8 x(r) 9 x(r) 10 x(r) 11 x(r) 12 x(r) 13 0.049 0.006 0.013 0.217 0 0 0 0 0.012 0.001 0.938 0.716 0.742 0.997 0.046 0.061 0.069 0.37 0.065 0.048

BCBM 0.613 0.024 0.036 0.915 0.006 0.003 0.005 0.053 0.002 0.001 Variable x(m) 1 x(m) 2 x(m) 3 x(m) 4 x(m) 5 x(m) 6 x(m) 7 x(r) 1 x(r) 2 x(r) 3 forward 0.219 0.216 0.258 0.954 0.208 0.204

1 0.011 0.013 0.012 0.955 0.945 0.991 0.996 0.993 0.995

1 0.462 0.504 0.642

From table 1 it can be seen that the forward selection often does not include the redundant predictors. However, the final subset almost does not get a significant part of the relevant features. Only x(m) , x7(m) correlated with the other significant predictors are constantly included in the subset.

4 This has a negative effect on the objective function values, which are far from optimal. The evaluation function values are displayed as a boxplot in figure 1. From figure 1, it can be concluded that the variation of the objective function values in the simulations is very large, and the average value is much smaller than that achieved by the backward elimination and the binary cut-and-branch method.

Both the backward elimination and the binary cut-and-branch method include all relevant features in the subset of predictors. But the backward elimination often leaves too many redundant predictors in the subset. Compared to this, for the BCBM results, such cases are relatively rare. Only a few redundant attributes x4(r) , x(r) are often included in the subset. This is because they are related to the 7 relevant predictors x(m) , x7(m) which are correlated with the other significant predictors. Nevertheless, 4 despite this problem, the binary binary cut-and-branch method provides the best values of the objective function compared to the backward elimination (figure 1).

The multicollinearity problem described is evidently difficult for the CFS method. It is a problem of the method itself, rather than optimization algorithms. This may have a negative effect on its applicability in practice.

5. Selection of leading indicators of consumer price index growth

Let us verify the applicability of the Correlation-based Feature Selection method using a real example. For this purpose, the task of forecasting the consumer price index was chosen. From economic research [14] it is known that one of the leading indicators of price changes during a business cycle is the industrial materials price index. The data provided by the unified interdepartmental informationalstatistical system [15] include the monthly time series of the price index for the acquisition of machines and equipment for investment purposes. Indices are grouped by type of economic activity and regions of the Russian Federation. Further, the territory of the Russian Federation as a whole is selected. 2 1 1 0 1 1 6 0 1 4 0 1 2 0 1 real ARIMA naive

The period from June to November 2010 was chosen as the test time interval. During that period, there was a sharp increase in the base price index in percent compared to the corresponding period of the previous year. The training dataset was taken from January 2006 to May 2010. It was used to fit the ARIMA model. But it naturally predicts the continuing fall in prices, as shown in figure 2. Let us verify whether this forecast can be improved by using leading indicators.

For the time period used to train the model, there are the data on price indices for the acquisition of machines and equipment for investment purposes for 92 types of economic activity, including 16 sections of Classification and a total index for all types of activity. The values of indices are in percent compared to the corresponding period of the previous year.

First, a naive approach was used, suggesting that the best predictor is the total index for all activities. The CFS method was used for optimal lags selection only. Lags from 0 to -6 were considered as candidates. As a result, lags -4 and -6 are chosen. This approach led to a rather low value of the objective function F* (table 2).

Next, the price indices for the acquisition of machines and equipment for investment purposes by economic activity, taken with lags from 0 to -6, were considered as candidate predictors for the application of greedy algorithms. Thus, the total number of predictors n was 644. The obtained optimal number n* is presented in table 2.

Method naive forward selection backward elimination forward selection backward elimination binary cut-and-branch

Finally, the number of predictors was reduced to indices by sections and the total index, that is, 17 indices. The total number of variables was 119 taking into account the possible lags. Table 2 shows the optimal values of the objective functions F* and the number of predictors n* obtained using the methods of forward selection, backward elimination, and binary cut-and-branch.

In both cases, the forward selection method leaves very few predictors. The values of the objective function are the lowest. The method of backward elimination selects a lot of variables, which, of course, are redundant. The binary cut-and-branch method provides the best value of the objective function.

In order to compare the forecasting performance on test data, regression models were estimated on the training set with the inclusion of selected variables with selected lags. For this, the dynlm package of the statistical environment R was used.

The future predictor values should be available in order to build a forecast for six months ahead. It was assumed that their real values are available only until May 2010. If later values were needed, then the predictor values were forecasted based on the ARIMA model. For the automatic selection of the structure of the ARIMA model, the auto.arima function of the forecast package was used. It is implemented in the R environment. For forecasting, the forecast function from the same package was used. Forecasting based on the results of the predictor selection by the backward elimination method was not performed, since the number of variables (39 and 18) is clearly redundant for estimating a model of 47 months (taking into account the earliest lag -6).

Table 3 presents the deviations of the real values of the consumer price index from its forecasts. The smallest absolute differences are in bold. It is revealed that the naive approach gives a good result. Graphically, it is shown in figure 2.

Method ARIMA naive forward selection, 92 indices forward selection, 17 indices binary cut-and-branch, 17 indices

This means that the Correlation-based Feature Selection approach can be recommended for choosing the optimal lags of predictor variables in a time series model. The simultaneous selection of price indices for the acquisition of machines and equipment for investment purposes by economic activity and their lags is more complicated. The time series of indices are very similar (some are even almost identical). Hence there is a very high correlation between the variables. At the same time, the lagged values of the indices are also highly correlated. This is similar to the case of multicollinearity from the model example considered in the previous section. As it was revealed above, in such a model, the achievement of optimum in problem ( 1 ) does not guarantee that only relevant features will be selected. In the structure of the solution, a certain proportion of redundant predictors is allowed.

Evidently, this problem also occurs in the selection of indices as leading indicators. As a result, the optimal solution obtained using the binary cut-and-branch method provides the best forecast performance for the medium term (2-3 months) only. In the long term, 5-6 months ahead, forecasts are worse than when selecting optimal lags for the total index only.

This effect is clearly visible when compared the results of selection by forward selection and binary cut-and-branch method for 17 indices. The results of forward selection give a smaller number of predictors and better predict for the long term. This can be explained by the overfitting effect, since using the binary cut-and-branch method, a greater number of predictors have been selected, some of which may be redundant.

6. Conclusions and recommendations

Thus, the Correlation-based Feature Selection method is applicable for selecting leading indicators in the time series forecasting. The candidate predictors are highly correlated. From simulation studies, it is revealed that greedy heuristics optimization algorithms in such cases do not give satisfactory results. The forward selection does not include many relevant features; the backward elimination leaves many redundant predictors in the subset. The binary cut-and-branch method gives the best result. However, the CFS method is not perfect: the optimal value of the heuristics does not guarantee that only all relevant predictors will be selected, and the inclusion of redundant variables is possible. This happens when the relevant features correlate with each other and with redundant variables. To avoid this, it is recommended by forming the initial set of candidate predictors to pre-exclude duplicate indicators and indicators with very similar dynamics.

Acknowledgments

The reported study was funded by Russian Ministry of Education and Science, according to the research project No. 2.2327.2017/4.6.

References

[1] Gaidel A V and Krasheninnikov V R 2016 Feature selection for diagnosing the osteoporosis by femoral neck X-ray images Computer Optics 40 ( 6 ) 939 - 946 DOI: 10.18287/ 2412 -6179-2016- 40-6- 939 -946

[2] Kutikova V V and Gaidel A V 2015 Study of informative feature selection approaches for the texture image recognition problem using the Laws' masks Computer Optics 39 ( 5 ) 744 - 750 DOI: 10.18287/ 0134 -2452-2015-39-5- 744 -750

[3] Bolón-Canedo V , Sánchez-Maroño N and Alonso-Betanzos A 2015 Knowledge-Based Systems 86 33 - 45

[4] Sagaert Y R , Aghezzaf E H , Kourentzes N and Desmet B 2018 European Journal of Operational Research 264 558 - 569

[5] Tibshirani R 1996 Journal of the Royal Statistical Society. Series B (Methodological) 58 267 - 288

[6] Fite J T , Don Taylor G, Usher J S , English J R and Roberts J N 2002 International Journal of Physical Distribution & Logistics Management 32 299 - 308

[7] Flach P 2012 Machine learning: the art and science of algorithms that make sense of data (Cambridge University Press)

[8] Hall M A 1999 Correlation-based feature selection for machine learning PhD thesis (Hamilton: University of Waikato)

[9] Nguyen H , Franke K and Petrovic S 2009 Proc. of the NIPS 2009 Workshop on Discrete Optimization in Machine Learning: Sub modularity, Sparsity & Polyhedra (DISCML) (Vancouver, Canada)

[10] Chang C-T 2001 European Journal of Operational Research 131 224 - 227

[11] Mezentsev Y A 2016 Binary Cut-and-Branch Method for Solving Linear Programming Problems with Boolean Variables CEUR Workshop Proceedings 1623 72 - 85

[12] Mezentsev Y 2017 Constructive Nonsmooth Analysis and Related Topics (dedicated to the memory of V.F . Demyanov ) (CNSA) (St. Petersburg , Russia) 1 - 3

[13] Sutter J M and Kalivas J H 1993 Microchemical journal 47 60 - 66

[14] Klein P A and Moore G H 1983 Journal of forecasting 2 119 - 135

[15] Unified interdepartmental informational-statistical system URL: https://fedstat .ru/