<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/combining-and-comparing-models-using-model-confidence-set/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://insightr.wordpress.com/2017/09/03/combining-and-comparing-models-using-model-confidence-set/"> R – insightR</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<div><strong>By Gabriel Vasconcelos</strong></div>
</div>
<div></div>
<div></div>
<p>In many cases, especially if you are dealing with forecasting models, it is natural to use a large set of models to forecast the same variable and then select the best model using some error measure. For example, you can break the sample into a training sample (in-sample) and a test sample (out-of-sample), estimate all models in the training sample and see how the perform in the test sample. You could compare the models using the root mean squared error (RMSE) or the mean absolute error (MAE).</p>
<p><span></span></p>
<p>If you want to make things a little more formal, then you can use some statistical test such as Diebold-Mariano to compare the models pairwise. This test will tell you if model A is better than model B. If you have <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> models to compare you will need to perform <img alt="n(n-1)/2" title="n(n-1)/2" src="https://s0.wp.com/latex.php?latex=n%28n-1%29%2F2&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> tests. If  <img alt="n" title="n" src="https://s0.wp.com/latex.php?latex=n&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is large you will have a lot of test results to look at and reporting the results will become a problem (imagine that you have 100 models, you will have 4950 p-values to report).</p>
<p>Fortunately, Hansen, Lunde and Nason proposed (<a rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.3982/ECTA5771/abstract" target="_blank">here</a>) a suitable solution in an article in 2011. Their solution is called Model Confidence Set (MCS). The MCS has an interpretation similar to a confidence interval, i. e. the <img alt="MCS_{1-alpha}" title="MCS_{1-alpha}" src="https://s0.wp.com/latex.php?latex=MCS_%7B1-%5Calpha%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> contains the best model with <img alt="100(1-alpha)%" title="100(1-alpha)%" src="https://s0.wp.com/latex.php?latex=100%281-%5Calpha%29%5C%25&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> confidence. The number of models in the MCS increases as we decrease <img alt="alpha" title="alpha" src="https://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> just like the size of a confidence interval. You only need to perform the MCS procedure once to compare all models and you have to report only the p-value for each model, which is an output of the test.</p>
<h2>Application</h2>
<p>Before going to the application, let us define a function that calculates the RMSE for several models at the same time, which we will be using a lot.</p>
<pre>rmse=function(real,...){
  forecasts=cbind(...)
  sqrt(colMeans((real-forecasts)^2))
}
</pre>
<p>To the application!! We are going to generate some data with the following design:</p>
<ul>
<li>300 observations (<img alt="t=1,dots,300" title="t=1,dots,300" src="https://s0.wp.com/latex.php?latex=t%3D1%2C%5Cdots%2C300&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">) and 8 variables (<img alt="n=1,dots,8" title="n=1,dots,8" src="https://s0.wp.com/latex.php?latex=n%3D1%2C%5Cdots%2C8&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">),</li>
<li>All variables are generated from <img alt="y_{i,t}=phi y_{i,t-1}+lambda_if_t +varepsilon_t" title="y_{i,t}=phi y_{i,t-1}+lambda_if_t +varepsilon_t" src="https://s0.wp.com/latex.php?latex=y_%7Bi%2Ct%7D%3D%5Cphi+y_%7Bi%2Ct-1%7D%2B%5Clambda_if_t+%2B%5Cvarepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> with <img alt="phi=0.6" title="phi=0.6" src="https://s0.wp.com/latex.php?latex=%5Cphi%3D0.6&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">, <img alt="lambda" title="lambda" src="https://s0.wp.com/latex.php?latex=%5Clambda&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is sampled from a uniform(-1,1) distribution and <img alt="varepsilon_{i,t}" title="varepsilon_{i,t}" src="https://s0.wp.com/latex.php?latex=%5Cvarepsilon_%7Bi%2Ct%7D&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is Normal with mean equals 0 and variance equals 4,</li>
<li><img alt="f_t" title="f_t" src="https://s0.wp.com/latex.php?latex=f_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is a common factor generated from <img alt="f_{t}=rho f_{t-1}+epsilon_t" title="f_{t}=rho f_{t-1}+epsilon_t" src="https://s0.wp.com/latex.php?latex=f_%7Bt%7D%3D%5Crho+f_%7Bt-1%7D%2B%5Cepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> with <img alt="rho=0.75" title="rho=0.75" src="https://s0.wp.com/latex.php?latex=%5Crho%3D0.75&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> and <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex"> is Normal with mean equals 0 and variance equals 4,</li>
<li>I generated two independent factors: four variables are generated from factor 1 and the remaining four variables are generated from factor 2.</li>
<li>We assume that we do not know which factor generates each variable.</li>
</ul>
<pre>library(reshape2)
library(grid)
library(gridExtra)
library(ggplot2)
library(boot)
T = 300 # = number of observations = #
n = 8 # = Number of Variables = #
rho = 0.75 # = factors AR coefficient = #
phi = 0.6 # = Variables AR coefficient = #

set.seed(1) # = seed for replication = #

# = Generate two autoregressive factors = #
f = matrix(0, T, 2)
for(i in 2:T){
  f[i, ] = rho * f[i-1, ] + rnorm(2)
}

lambda = runif(n,-1,1) # = Factor loadings = #

# = Generate variables = #
Y = matrix(0, T, n)
for(i in 2:T){
  # = Based on the first factor = #
  Y[i, 1:4] = phi * Y[i-1, 1:4] + lambda[1:4] * f[i, 1] + rnorm(4, 0, 2)
  # = Based on the second factor = #
  Y[i, 5:8] = phi * Y[i-1, 5:8] + lambda[5:8] * f[i,2] + rnorm(4, 0, 2)
}

# = Organize variable lags = #
# = y is the first (dependent) variable = #
Y = as.data.frame(embed(Y, 2)[ , c(1, 9:16)])
# = y is the first (dependent) variable = #
colnames(Y) = c(&quot;y&quot;, paste(&quot;x&quot;, 1:n, sep = &quot;&quot;))
</pre>
<p>Now that we have our data, we are going to estimate all possible combinations of linear regressions using our 8 variables. The first variable will be the dependent variable, defined as <img alt="y" title="y" src="https://s0.wp.com/latex.php?latex=y&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. The explanatory variables are the first lags of all eight variables (including the first), defined as <img alt="x_i, i=1,dots,8" title="x_i, i=1,dots,8" src="https://s0.wp.com/latex.php?latex=x_i%2C+i%3D1%2C%5Cdots%2C8&bg=ffffff&%23038;fg=333333&%23038;s=0" class="latex">. This design results in 255 models to be estimated and compared.</p>
<pre># = First we will generate formulas for all the 255 equations = #
variables = colnames(Y)[-1]
combinations = lapply(1:n, function(x)combn(1:n, x))
formulas = lapply(combinations, function(y){
  apply(y, 2, function(x){
    paste(&quot;~&quot;, paste(variables[x], collapse = &quot; + &quot;), sep=&quot; &quot;)
  })
})
formulas = unlist(formulas)

# = Now we break the sample into in-sample and out-of-sample = #
in_sample = Y[1:100, ]
out_of_sample = Y[-c(1:100), ]

# = Estimate the models = #
models = lapply(formulas, function(x) lm(paste(&quot;y&quot;, x), data = in_sample))
# = Compute forecasts = #
forecasts = Reduce(cbind, lapply(models, function(x) predict(x, out_of_sample)))
colnames(forecasts) = paste(&quot;model&quot;, 1:ncol(forecasts), sep = &quot;&quot;)

# RMSES (Mean of all models, model with all variables)
rmse(out_of_sample$y,rowMeans(forecasts),forecasts[,ncol(forecasts)])
</pre>
<pre>## [1] 2.225175 2.322508
</pre>
<p>As you can see, the average forecast is more accurate than the model using all 8 variables. Now we are going to compute the MCS using the function below. If you want more details on the algorithm, just look at the original article <a rel="nofollow" href="http://onlinelibrary.wiley.com/doi/10.3982/ECTA5771/abstract" target="_blank">here</a>.</p>
<pre>mcs=function(Loss,R,l){
 LbarB=tsboot(Loss,colMeans,R=R,sim=&quot;fixed&quot;,l=l)$t
 Lbar=colMeans(Loss)
 zeta.Bi=t(t(LbarB)-Lbar)

  save.res=c()
  for(j in 1:(ncol(Loss)-1)){
    Lbardot=mean(Lbar)
    zetadot=rowMeans(zeta.Bi)
    vard=colMeans((zeta.Bi-zetadot)^2)

    t.stat=(Lbar-Lbardot)/sqrt(vard)
    t.max=max(t.stat)
    model.t.max=which(t.stat==t.max)

    t.stat.b=t(t(zeta.Bi-zetadot)/sqrt(vard))
    t.max.b=apply(t.stat.b,1,max)

    p=length(which(t.max&lt;t.max.b))/R

    save.res=c(save.res,p)
    names(save.res)[j]=names(model.t.max)

    Lbar=Lbar[-model.t.max]
    zeta.Bi=zeta.Bi[,-model.t.max]

  }
  save.res=c(save.res,1)
  names(save.res)[j+1]=names(Lbar)
  save.p=save.res
  for(i in 2:(length(save.res)-1)){
    save.p[i]=max(save.res[i-1],save.res[i])
  }
  aux=match(colnames(Loss),names(save.p))
  save.p=save.p[aux]
  save.res=save.res[aux]
  return(list(test=save.p,individual=save.res))
}
</pre>
<p>The MCS is computed in a single line of code, but first we are going to break the out-of-sample into two sub-samples. We are going to compute the MCS in the first sub-sample and compare the models in the second sub-sample. The MCS function returns the p-value needed to include each model in the set. For example, if model1 has a p-value=0.1, it will be included in the MCS of 90\% confidence.</p>
<pre># = generate 2 subsamples = #
out1 = out_of_sample[1:100, ]
out2 = out_of_sample[-c(1:100), ]
f1 = forecasts[1:100, ]
f2 = forecasts[-c(1:100),]

# = Compute MCS = #
set.seed(123) # = Seed for replication = #
MCS=mcs( Loss=(f1-out1$y)^2, R=1000, l=3 )$test # = R block bootstrap samples with block length = 3  = #
# = If you want to see the p-values write (print(sort(MCS)) = #
# = Selected models in the 50% MCS  = #
selected=which(MCS>0.5)

# RMSE (Mean of all models, model with all variables, MCS models)
rmse(out2$y,rowMeans(f2),f2[,ncol(f2)],rowMeans(f2[,selected]))
</pre>
<pre>## [1] 2.142932 2.240110 2.013547
</pre>
<p>As you can see, if we combine only the models in the MCS the RMSE becomes smaller. This is not always the case, but it happens a lot. The combined results are better when the forecasts do not have a very high positive correlation. The MCS is not useful only to combine forecasts. Even if the combinations are not better than individual models, you still have a confidence set that includes the best model with the desired significance. If some model is excluded from the confidence set you have strong evidence that this model is worst than the ones in the set (depending on the confidence level, of course). Naturally, if you increase the confidence level, the number of models in the set also increases just like in confidence intervals.</p>
<p>Finally, we can find some interesting relationships in the plots below. The first plot shows the forecasting errors as function of the confidence level. The optimal confidence level is around 55\%. If we increase it more the MCS starts to select a lot o bad models and the forecast becomes less accurate. The second plot shows the number of models as a function of the confidence level and the third plot shows the errors as a function of the number of models.</p>
<pre>alpha = seq(0,1,0.01)
nmodels = sapply(alpha, function(x) length(which(MCS >= x)))
errors = sapply(alpha, function(x) rmse(out2$y, rowMeans(as.matrix(f2[ , which(MCS >= x)]))))

df = data.frame(confidence = 1 - alpha, errors = errors, nmodels = nmodels)
dfmelt = melt(df, id.vars = &quot;errors&quot;)

g1=ggplot(data=df) + geom_point(aes(x=confidence,y=errors))
g2=ggplot(data=df) + geom_point(aes(x=confidence,y=nmodels))
g3=ggplot(data=df) + geom_point(aes(x=nmodels,y=errors))
grid.arrange(g1,g2,g3,ncol=2)
</pre>
<p><img alt="plot of chunk unnamed-chunk-6" src="https://i1.wp.com/i.imgur.com/3hSwHd8.png?w=456&ssl=1" data-recalc-dims="1"></p>
<h2>References</h2>
<p>Hansen, Peter R., Asger Lunde, and James M. Nason. “The model confidence set.” Econometrica 79.2 (2011): 453-497.</p>
<p> <a rel="nofollow" href="http://feeds.wordpress.com/1.0/gocomments/insightr.wordpress.com/650/" target="_blank"><img border="0" alt="" src="http://feeds.wordpress.com/1.0/comments/insightr.wordpress.com/650/"></a> <img border="0" alt="" src="https://i1.wp.com/pixel.wp.com/b.gif?resize=1%2C1&ssl=1" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://insightr.wordpress.com/2017/09/03/combining-and-comparing-models-using-model-confidence-set/"> R – insightR</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/7sZylKl8AOo" height="1"></a></div>