Forecasting with Model Uncertainty: Representations and Risk
Reduction*
Keisuke Hirano†

Jonathan H. Wright‡

First version: October 14, 2013This version: April 23, 2016

Abstract
We consider forecasting with uncertainty about the choice of predictor variables. The researcher
wants to select a model, estimate the parameters, and use the parameter estimates for forecasting.
We investigate the distributional properties of a number of different schemes for model choice andparameter estimation: in-sample model selection using the Akaike information criterion, out-ofsample model selection, and splitting the data into subsamples for model selection and parameter
estimation. Using a weak-predictor local asymptotic scheme, we provide a representation result thatfacilitates comparison of the distributional properties of the procedures and their associated forecast risks. We develop a simulation procedure that improves the accuracy of the out-of-sample andsplit-sample methods uniformly over the local parameter space. We also examine how bootstrap aggregation (bagging) affects the local asymptotic risk of the estimators and their associated forecasts.
Numerically, we find that for many values of the local parameter, the out-of-sample and split-sample
schemes perform poorly if implemented in the conventional way. But they perform well, if implemented in conjunction with our risk-reduction method or bagging.

* We are grateful to Don Andrews, Marine Carrasco, Russell Davidson, Gary Chamberlain, Sylvia Gonçalves, Bruce Hansen,
Serena Ng, Peter Phillips, Jack Porter, James Stock, Mark Watson, four anonymous referees and the co-editor for very helpful
discussions. The usual disclaimer applies.
† Department of Economics, University of Arizona, 1130 E.Helen St., Tucson AZ 85721. Email: hirano@u.arizona.edu‡ Department of Economics, Johns Hopkins University, 3400 North Charles St., Baltimore MD 21218.
Email:

wrightj@jhu.edu

11 Introduction
In this paper, we reconsider the problem of forecasting when there is uncertainty about the forecastingmodel. It is common practice to select the model based on pseudo-out-of-sample fit from a sequence ofrecursive or rolling predictions. Parameters are then estimated over the whole sample period. The ideaof using an out-of-sample criterion is discussed by Clark (2004) and West (2006), and has a long history,
going back to Wilson (1934), Meese and Rogoff (1983), and Ashley, Granger, and Schmalensee (1980).
Instead, one might select the model based on in-sample fit, but adjust for potential overfitting by using
an information criterion, such as the Akaike Information Criterion (AIC) (Akaike, 1974), as advocated byInoue and Kilian (2006).
We consider a pseudo-likelihood setting with a fixed number k of potential parameters to be estimated,each of which has a coefficient that is local to zero. The concept of model selection that we envision
amounts to selecting a set of zero restrictions; in a regression setting, for example, this would indicate
which predictors are excluded from the regression. Thus there are up to 2k possible models among which
we can choose. Having chosen the model, we then have to estimate the parameters and use these forforecasting. Although some model will be best in terms of predictive accuracy, the local-to-zero nestingmeans that we can never consistently select that model. We consider various methods of model selection
and forecasting, including: using in-sample fit with the AIC information criterion; selecting the model
based on recursive pseudo-out-of-sample forecast accuracy and then using the whole dataset for parameter estimation; and splitting the sample into two parts, using one part for model selection and theother for parameter estimation. We call this last method the split-sample approach. Unlike the first two
methods, it is not commonly used in practice. But it does ensure asymptotic independence between
parameter estimates and model selection, unlike methods based on in-sample fit (Leeb and Pötscher,
2005; Hansen, 2009), and also unlike the standard out-of-sample approach.
We obtain asymptotic characterizations of these forecasting procedures under the local parameter sequence. A key step is to obtain an asymptotic representation of the partial sum process for the scorefunction as the sum of a term that is directly informative about the local parameters, and another termthat is an independent Gaussian process. This allows us to provide a limit-experiment type representation of the procedures, from which we can calculate normalized local asymptotic mean square prediction

2

errors up to O(T −1 ) terms. We show that the recursive pseudo-out-of-sample and split-sample procedures are inefficient, in the sense that their limit distributions depend on the ancillary Gaussian noise
process. Our characterizations also suggest ways to improve upon these procedures. The influence of
the ancillary noise term in the limiting asymptotic representation can be eliminated by a conditioningargument. We can implement this noise reduction via a simulation-and-averaging scheme; doing thisis shown to uniformly improve the out-of-sample and split-sample methods asymptotically for a wide
variety of loss functions.
This method is related to bootstrap aggregating (bagging) (Breiman, 1996) in which the data are resampled, the forecasting method is applied to the resampled data, and the resulting forecasts are then averaged over all the bootstrap samples. Bagging has a smoothing effect that alters the risk properties of estimators, but averaging over bootstrap draws can also reduce the influence of the extraneous noise termin the out-of-sample and split-sample methods. Earlier theoretical work on bagging, notably Bühlmann
and Yu (2002), emphasized its smoothing effect rather that the noise reduction effect, though Efron(2014) observes that it can reduce variability. 1
We then numerically compare the various procedures in terms of their local asymptotic risk. In theirstandard forms the rank ordering among the in-sample, out-of-sample and split-sample methods depends on the localization parameter which is not consistently estimable. Nonetheless, we find that for
many values of the localization parameter, in-sample forecasting using the AIC gives the most accurate
forecasts, out-of-sample prediction does worse, and the split-sample method does worst of all. This
is intuitive because the out-of-sample and split-sample schemes are in some sense wasting data, andis essentially the argument of Inoue and Kilian (2004) and Inoue and Kilian (2006) for the use of insample rather than out-of-sample predictability tests. However, introducing the simulation or bagging
step changes the rank ordering substantially, and entirely reverses it for many values of the localization
parameter. Our simulation scheme has no effect on in-sample forecasts, but reduces the local asymptotic mean square prediction error of the out-of-sample and split-sample forecasts uniformly in the localization parameter, and the reductions are generally numerically large. In short, the out-of-sample andsplit-sample forecasts are poor choices in their conventional form, but become very competitive when
1 One other useful feature of the out-of-sample forecasting setup is that it can be constructed to use only real-time data which
precisely mimics the data available to a researcher in the presence of data revisions. Unfortunately, adding our simulation
scheme or bootstrap aggregation step destroys this feature.

3

combined with a simulation or bagging step.
In the next section, we set up the model and introduce the various procedures we will evaluate. In Section
3, we derive asymptotic characterizations via our representation theorem. Section 4 contains some extensions. Section 5 explores the local asymptotic risk properties of the procedures numerically. Section
6 contains an illustrative application, and Section 7 concludes.

2 Pseudo-Likelihood Framework
We observe (y t , x t ) for t = 1, . . . , T and wish to forecast y T +1 given knowledge of x T +1 . Let the pseudo log(conditional) likelihood be`(β) =TX`t (β) =

t =1

T
X

ln f (y t |x t , β),

(2.1)

t =1

where f is a conditional density function and the parameter β is k × 1. This framework could apply
to cross-sectional regression of an outcome variable y t on a k × 1 vector of predictors x t , h-step ahead
forecasting regressions (where x t are suitably lagged predictor variables), vector autoregression (wherex t contains lagged values of the vector y t ), and nonlinear regression models. There may be unmodeled
dependence, subject to the large-sample distributional assumptions imposed below.
We consider a setting where the models under consideration are nested in such a way that model selection amounts to setting some elements of β to zero, and estimating the others. Thus there are up to 2k
possible models. Let m ⊂ {1, . . . , k} denote a model, with the interpretation that the elements of m indicate which coefficients of β are allowed to be nonzero. The set of possible models M is a subset of thepower set of {1, . . . , k}.
We consider a variety of strategies for model selection and parameter estimation. Each strategy will end
up giving us an estimator β̃ for β, some elements of which may be zero. The strategies that we considerare:

1. Maximum Likelihood (MLE). Set β̃ to the unrestricted (pseudo) maximum likelihood estimator,β̂, that maximizes `(β).
2. James-Stein (JS). The positive-part James-Stein estimator uses the unrestricted estimate β̂ and an4

estimate V̂ of its asymptotic variance-covariance matrix. The JS estimator for k > 2 is:(β̃ = β̂ · max 1 −k −2
T β̂0V̂ −1 β̂

)

,0 .

(2.2)

3. AIC (In-sample). For each model m ∈ M , let β̂(m) denote the restricted pseudo-ML estimatorthat maximizes `(β) subject to the zero restrictions implied by m. (Thus, β̂(m) is a k × 1 vector
with zeros in the locations implied by m.) Let n(m) be the number of free parameters in model m.
For each m ∈ M we calculate the AIC objective function

AIC (m) = 2`(β̂(m)) − 2n(m),(2.3)

and choose the model m ∗ that maximizes AIC (m). Then set β̃ = β̂(m ∗ ).
4. Out of Sample. For each model m, we estimate the model recursively starting a fraction π ∈ (0, 1)of the way through the sample, and calculate its one-period-ahead predictive density to obtain apseudo out-of-sample estimate of predictive performance. Let β̂1,t −1 (m) denote the pseudo maximum likelihood estimate for model m using observations 1 to t − 1. For each m, we calculateTX
t =[T π]+1`t (β̂1,t −1 (m)).

(2.4)

We then choose the model m that maximizes this predictive likelihood, and use the full sample for
estimation of the model.2
5. Split-sample. For each model m, we calculate AIC using data up to a fraction π of the way throughthe sample:
AIC ss (m) = 2`(β̂1,[T π] (m)) − 2n(m).

(2.5)

For m ∗ = arg max AIC ss (m), we use the second fraction of the sample to estimate the model parameters:β̃ = β̂[T π]+1,T (m ∗ ).

(2.6)

2 Several authors test for the statistical significance of differences in out-of-sample forecasting performance with one model,

typically a simple benchmark, as the null (see, for example Diebold and Mariano (1995) and Hansen and Timmermann (2015)).
Here we are instead thinking of selecting the model based on the point estimate of its pseudo out-of-sample predictive performance.

5Later in the paper, we also consider modifications of these procedures that attempt to improve their risk
properties.
We focus on obtaining limiting distributions for these estimators, which we denote generically as β̃. Thismay be of direct interest in itself as a parameter estimation problem. Alternatively, it may be usefulfor forecasting. Suppose that a forecaster creates a point forecast of y T +1 using a rule ŷ T +1 (β̃). The
function ŷ T +1 (·) may depend on x T +1 , and it may depend on the data up to time T in other ways. Given
a forecasting loss function L( ŷ T +1 , y T +1 ), let£¤L(a, β) = E T L( ŷ T +1 (a), β) ,(2.7)

where E T [.] denotes the expectation conditional on the data up to and including time T . Define theregret as:
L r (a, β) = L(a, β) − inf L(c, β).

(2.8)

c

This subtracts off the loss associated with the infeasible optimal choice of the estimate used in the forecast rule. The regret risk is defined as:hi
R(β̃, β) = E L r (β̃, β) ,(2.9)

which depends on the “true” parameter value β and the estimator β̃ (the rule for estimation, not the
sample value).
Typically, our estimators β̃ will bep
T -consistent for β, so a suitably normalized version of the forecast

risk will depend on the local behavior of L r (a, β) as a → β. For a number of standard forecast loss functions, the normalized regret risk for a suitably chosen forecast rule ŷ T +1 (β̃) will satisfy:h ³p´i
T · R(β̃, β) = E G T (β̃ − β) + o (1) .

for some convex function G(·). Thus obtaining the limiting distribution of(2.10)p
T (β̃ − β) will be key to char-acterizing the risk properties of the associated forecast rule. The fact that the regret risk is convex inp
T (β̃ − β) will also prove to be useful for our risk reduction methods. In the next section we give someexamples of some cases where the representation in Equation (2.10) holds.
62.1 Example: Prediction in a Regression Model
To illustrate our approach in a simple setting, we consider prediction using a standard regression model:

y t = β0 x t + u t ,(2.11)

where the u t are i.i.d. with mean 0, finite variance σ2 , and 2+δ finite moments for some δ > 0. We assume
that x t is a k × 1 i.i.d. vector that has been orthonormalized, so that E [x t x t0 ] = I k .3
This model fits into the general pseudo-likelihood framework of Section 2, using the standard Gaussian
likelihood. Then β̂, the unrestricted pseudo-ML estimator of β, is the OLS estimator; and β̂(m), therestricted pseudo-ML estimator under model m, is the restricted OLS estimator using only the regressors
indicated by m.
Consider forecasting under squared error loss, with forecasts of the form ŷ T +1 = β̃0 x T +1 , where β̃ canagain be any of the estimators of β. The mean squared prediction error is£¤£¤
M SPE = E (y T +1 − β̃0 x T +1 )2 = E (u T +1 − (β̃ − β)0 x T +1 )2£¤= σ2 + E (β̃ − β)0 (β̃ − β) .

(2.12)

The first term on the right hand side of (2.12) is the risk of infeasible forecast rule that uses the true
value of the parameter, β0 x T +1 . The second term is O(T −1 ) and differs across forecasting methods. In
the forecasting context, parameter estimation uncertainty is of small order relative to uncertainty in the
distribution of the shocks. But uncertainty about the distribution of the shocks is common to all themodels. We therefore subtract the term σ2 and scale the risk by T to obtain the normalized mean squareprediction error criterion
i
hp
p
¡
¢
T (β̃ − β)0 T (β̃ − β) .
N M SPE = T M SPE − σ2 = E(2.13)

This is equivalent to the normalized regret risk, as defined in equation (2.10), for squared error loss.
Other loss functions lead to different expressions for the regret risk, but in some leading cases they will
3 The orthonormality of x is not essential for the analysis, but simplifies the notation.
t7

continue to satisfy equation (2.10). For example suppose we work with linex forecast loss,¡¢L( ŷ T +1 , y T +1 ) = exp θ(y T +1 − ŷ T +1 ) − θ(y T +1 − ŷ T +1 ) − 1,

(2.14)

for some given θ 6= 0, and suppose that u T +1 is Gaussian. In this case, the oracle point forecast using
knowledge of β and σ2 is y T∗ +1 = β0 x T +1 + θσ
2 , so it is natural to consider rules of the form
2ŷ T +1 (β̃) = β̃0 x T +1 +

θ σ̂2
,
2

where σ̂2 is a consistent estimator of σ2 . This forecast rule employs an estimate of the optimal “shift” ofthe naive forecast to account for the asymmetry of the linex loss function (Christoffersen and Diebold,1997). Using a second-order Taylor expansion, it can be shown that the normalized regret risk for thisforecast rule is
T · R(β̃, β) =ip
θ 2 hpE
T (β̃ − β)0 Var(x T +1 ) T (β̃ − β) + o(1),2(2.15)

again satisfying equation (2.10). With mean absolute loss, the normalized regret risk for forecasts of theform ŷ T +1 = β̃0 x T +1 also has the form in equation (2.10):

T · R(β̃, β) =ipf u (0) hp
E
T (β̃ − β)0 T (β̃ − β) + o(1),2(2.16)

where f u (·) is the probability density function of u T +1 .
To gain further intuition for our theoretical results in the next section, consider the special case whereP
the u t are i.i.d. N (0, 1) and the regressors are treated as fixed and satisfy T1 Tt=1 x t x t0 = I k . Then the¡¢least squares estimator for the full set of parameters β̂ is distributed N β, σ2 I k /n , and β̂ is a minimal
sufficient statistic for β. If a procedure makes nontrivial use of information in the data other than that
contained in β̂, it is introducing an unnecessary source of randomness. In the next section we will obtainan asymptotic analog to this argument in the general pseudo-likelihood setting, and show how it applies
to the various procedures we consider.

8

3 Local Asymptotics
In order to capture the role of parameter and model uncertainty in our analysis, the joint distribution of{(y 1 , x 1 ), . . . , (y T , x T )} is assumed to be a triangular array with drifting parameters. Let {(y 1 , x 1 ), . . . , (y T , x T )}have joint distribution P T , and define the pseudo-true value of the parameter as

β0,T = arg max
βZ`(β)d P T .

(3.1)

We assume that there is a unique solution to the maximization problem, and take the pseudo-true values(or functions of them) as our objects of interest. We will take limits as T → ∞ under the assumption thatbβ0,T = p ,Tb ∈ Rk .

(3.2)

This type of drifting sequence has been used by Claeskens and Hjort (2008) and Inoue and Kilian (2004)
and others to study the large-sample properties of model selection procedures. It preserves the role ofparameter uncertainty in the asymptotic approximations, unlike fixed-alternative asymptotics in whichmodel selection can determine which coefficients are nonzero with probability approaching one. Laterin the paper, we will show some Monte-Carlo evidence indicating that this local parameterization provides a good approximation in small samples. The analysis could be extended to allow some componentsof β to be localized away from zero, corresponding to situations where some components of β are known
to be nonzero. We use →d to denote weak convergence and →p to denote convergence in probability under the sequence of measures {P T }∞T =1 . Our results below depend crucially on the convergence propertiesof the partial sums of the pseudo-likelihood function. We make the following high level assumptions.

Assumption 3.1
T −1/2[Tr
X] ∂`t (β0,T )
t =1

∂β

→d B (r ),where B (r ) is a k-dimensional Brownian motion with positive definite covariance matrix Λ.
Assumption 3.2 For all sequences βT in a T −1/2 -neighborhood of zero,

−T −1

[Tr
X] ∂2 `t (βT )→p r Σ.
0
t =1 ∂β∂β

9where Σ is positive definite.

Assumption 3.3 For any fixed k × k matrix C ,0[Tr−1 ∂l (βX] tX
s 0,T )
t =1 s=1∂β∂l t (β0,T )
1
C
→d
t −1
∂β

Z r
0

1
B (s)C d B (s).
sThese high-level assumptions would follow from conventional regularity conditions in correctly specified parametric models. Assumptions 3.1 and 3.2 are time-series versions of standard results for thescore and Fisher information, while Assumption 3.3 is a technical condition that can often be shown,for example using results in Davidson (1994). In misspecified models, the assumptions require that thepseudo-true parameter sequence β0,T is related to the distribution of the data in a smooth way.
To gain intuition for the results to follow, consider the case where the parametric model with conditional
likelihood f (y t |x t , β) is correctly specified. Then, under standard regularity conditions, Assumptions 3.1,
3.2 and 3.3 will hold with Λ = Σ. Furthermore, the model will be locally asymptotically normal (LAN), andpossess a limit experiment representation (see for example van der Vaart 1998, Chs. 7-9). In particular,consider any estimator sequence β̃ with limiting distributions in the sense thatT 1/2 β̃ →d L b ,(3.3)

where the limit is taken under the drifting sequences of measures corresponding to β0,T = T −1/2 b, andL b is a law that may depend on b. Then the estimator β̃ has an asymptotic representation as a randomized estimator in a shifted normal model: if Y is a single draw from the N (Σb, Σ) distribution, and U
is random variable independent of Y (with sufficiently rich support4 ), there exists an estimator S(Y ,U )with
S(Y ,U ) ∼ L b(3.4)

for all b. In other words, the sequence T 1/2 β̃ is asymptotically equivalent to the randomized estimator S
under all values of the local parameter.
We extend this type of asymptotic representation, in terms of an asymptotically sufficient component
4 Typically, a representation S(Y ,U ) exists for U distributed uniform on [0, 1], but for our results below, it is useful to allow U
to have a more general form.

10and an independent randomization, to the pseudo-likelihood setup. We do this by establishing a largesample representation of the partial sum process for the score function that corresponds to the (Y ,U )limit experiments in parametric LAN models.
From Assumptions 3.1 and 3.2, it follows that:T −1/2[Tr
X] ∂`t (0)
t =1

∂β

→d B (r ) + r Σb =: Y (r )(3.5)

Thus the partial sums of the score function evaluated at β = 0 converge to a Brownian motion with linear
drift. By a standard argument, we can decompose this process into the sum of a normal random vector
and a Brownian bridge:

Proposition 3.4 Under the drifting sequence in Equation (3.2) and Assumptions 3.1 and 3.2,

T −1/2[Tr
X] ∂`t (0)
t =1

∂β

→d Y (r ) = r Y +UB (r ),(3.6)

where Y := Y (1) ∼ N (Σb, Λ), and UB (r ) is a k-dimensional Brownian bridge with covariance matrix Λ,
where UB is independent of Y .

All proofs are given in Appendix A. This result decomposes the limit of the partial sums of the scorefunction into two stochastic components, one of which depends on b and one of which is ancillary.
Let Σ(m) denote the k × k matrix that consists of the elements of Σ in the rows and columns indexed
by m and zeros in all other locations, and let H (m) denote the Moore-Penrose generalized inverse of
Σ(m). Then T 1/2 β̂ →d Σ−1 Y (1) and T 1/2 β̂(m) →d H (m)Y (1). This leads to the following asymptotic
characterizations of the procedures:

Proposition 3.5 Under the drifting sequence in Equation (3.2) and Assumptions 3.1, 3.2 and 3.3, we havethe following limiting representations of the parameter estimation procedures:(i) Using unrestricted MLE:
T 1/2 β̂ →d Σ−1 Y (1);

11

(3.7)

(ii) Using the positive-part James-Stein estimator:

T

1/2

½
β̃ →d Σ Y (1) max 1 −
−1

¾
k −2
,0 ;
Y (1)0 Σ−2 Y (1)(3.8)(iii) Selecting the model using the AIC:

T 1/2 β̃ →d

X

n
£
¤o
H (m ∗ )Y (1)1 m ∗ = arg max Y (1)0 H (m)Y (1) − 2n(m) ;

(3.9)

m

m∗

(iv) Selecting the model minimizing recursive out-of-sample error starting a fraction π of the way throughthe sample:

T 1/2 β̃ →dXnh
H (m ∗ )Y (1)1 m ∗ = arg max Y (1)0 H (m)Y (1)
m

m∗

io
1
− Y (π)0 H (m)Y (π) + t r (H (m)Λ) log(π) ;π(3.10)(v) Using the split-sample method, using the first fraction π of the sample for model selection and the restfor parameter estimation:

T

1/2

β̃ →d

X
m∗

where

∗ Y (1) − Y (π)

H (m )

1−π

½

·

1
1 m = arg max
Y (π)0 H (m)Y (π) − 2n(m)
m
π
∗

¸¾

;

(3.11)m ∗ denotes the summation over all the models in M .

P

Of course, there are other criteria besides AIC that we could use for in-sample model selection. Someof these are asymptotically equivalent to AIC, such as Mallows’ C p criterion (Mallows, 1973) or leaveone-out cross-validation. Using any of these information criteria for in-sample model selection will givethe same asymptotic distribution as in equation (3.9). Alternatively, one could use the Bayes information
criterion (BIC). In the present setting, because the penalty term goes to zero at a rate slower than T −1 , theBIC will pick the small model (β = 0) with probability converging to one. Part (v) of the proposition can
immediately be adapted to selecting the model minimizing out-of-sample error with a rolling estimation
window, as long as the estimation window contains a fixed fraction of the sample size, but not if it instead
contains a fixed number of observations as in Giacomini and White (2006).
Inoue and Kilian (2004) considered the local power of some in-sample and out-of-sample tests of the12hypothesis that β = 0. They derived equation (3.7) and a result very similar to equation (3.10). The
expression in equation (3.10) involves asymptotically selecting the model by maximizing an expressionthat is the difference between two quadratic forms plus a nonstochastic term. This derivation is also
closely related to Hansen and Timmermann (2015) who showed that the difference between the out-ofsample predictive accuracy of two models is asymptotically equivalent to the difference between two
Wald statistics plus a nonstochastic term.

3.1 Rao-Blackwellization
The estimators other than the out-of-sample and split-sample estimators can be viewed as generalized
shrinkage estimators (Stock and Watson, 2012) as their limiting distributions are of the form: T 1/2 β̃ →dY g (Y ) for some nonlinear function g (Y ). In contrast, the limiting distributions in equations (3.10) and(3.11) are functions of both Y and an independent Brownian bridge, UB (r ). Their dependence on thenoise term U = UB suggests a novel way to improve them.
In the statistical experiment of observing the pair (Y ,U ), where Y ∼ N (Σb, Λ) and U is ancillary, thevariable Y is sufficient. Thus, for any estimator S(Y ,U ), consider its conditional expectation

S̃(Y ) = E [S(Y ,U )|Y ].

(3.12)

By the Rao-Blackwell theorem, the risk of S̃(Y ) is less than that of S(Y ,U ) for all b for any strictly convex
risk function, including the normalized regret risk in equation (2.10), as long as S(Y ,U ) depends on U .
To implement the conditional estimators, we need consistent estimators Λ̂ →p Λ and Σ̂ →p Σ. Dependence in the scores poses no problem, so long as Λ̂ is a consistent estimate of the zero-frequency spectral
density. Take L independent artificially generated Brownian bridges {UBi (r )}Li=1 with covariance matrix
Λ̂. For each i , consider the estimators:

β̃i ,1 =Xnhβ̂(m ∗ )1 m ∗ = arg max T β̂(m)0 Σ̂β̂(m)
m

m∗

io
1
− [T 1/2 πβ̂(m) + Ĥ (m)UBi (π)]0 Σ̂[T 1/2 πβ̂(m) + Ĥ (m)UBi (π)] + t r (H (m)Λ̂) log(π)π

13(3.13)

and"β̃i ,2 =X∗β̂(m ) − T−1/2Ĥ (m)

Ubi (π)

#

1−π
½
·
¸¾
1
·1 m ∗ = arg max[πβ̂(m) + T −1/2 Ĥ (m)UBi (π)]0 Σ̂[πβ̂(m) + Ĥ (m)T −1/2UBi (π)] − 2n(m) ,m
πm∗(3.14)

where Ĥ (m) is the Moore-Penrose inverse of Σ̂(m). The next proposition gives their limiting distributions:

Proposition 3.6 Under the conditions for Proposition 3.5, for each i :T 1/2 β̃i ,1 →dXnh
H (m ∗ )Ỹi (1)1 m ∗ = arg max Ỹi (1)0 H (m)Ỹi (1)
m

m∗

io
1
− Ỹi (π)0 H (m)Ỹi (π) + t r (H (m)Λ) log(π)π(3.15)

andT 1/2 β̃i ,2 →dXm∗

H (m ∗ )½·¸¾
1
Ỹi (1) − Ỹi (π)
1 m ∗ = arg max
Ỹi (π)0 H (m)Ỹi (π) − 2n(m) ,m
1−ππ(3.16)

where Ỹi (r ) = r Y +UBi (r ). These are the same distributions as in equations (3.10) and (3.11).

These estimators can then be averaged over i . After this step of averaging over different realizationsof the Brownian bridge, the asymptotic distributions depend on Y alone and are asymptotically the
expectations of the out-of-sample and split-sample estimators conditional on Y . Note that this RaoBlackwellization (henceforth, RB) does applies only to the out-of-sample and split-sample estimators
because it is only for these estimators that there is any ancillary noise process to eliminate. But for these,
the limiting risk of the RB estimator is strictly lower than the original estimator, indicating that the outof-sample and split-sample estimators are asymptotically inadmissible.

3.2 Linear Regression Model and Bagging
In the special case of the regression model with orthonormal regressors, considered in Subsection 2.1,
we have Λ = Σ = σ−2 I k . In this model, all of the estimators depend crucially on the partial sum process14T −1/2

P[Tr ]

t =1 x t y t and it follows from Proposition 3.4 that:T −1/2 σ−2[TrX]x t y t →d Y (r ),(3.17)

t =1

and Proposition 3.5 will immediately apply.
In the linear regression model (subsection 2.1), we can also consider adding a bagging (bootstrap aggregation) step to each of the procedures (Breiman (1996)). Bühlmann and Yu (2002) studied the largesample properties of bagging. The i th bagging step resamples from the pairs {(x t , y t ), t = 1, . . . , T } withreplacement to form a pseudo-sample {x t∗ (i ), y t∗ (i ), t = 1, . . . , T }. The full model-selection and estimation procedure is then applied to the i th bootstrap sample. This is repeated L times, and the L estimates
are averaged to obtain the bagged estimate that can be used for forecasting. The following proposition
provides a key result for obtaining the limiting distribution with bagging.

Proposition 3.7 Let {x t∗ , y t∗ , t = 1, . . . , T } denote a bootstrap sample and g (·) be an uniformly integrable
Rk -valued functional. Then
" Ã

E

∗

X] ∗ ∗
1 [Tr
g
xt yt
p
σ2 T t =1!#´i
h ³ r p
→ E ∗ g 2 T β̂ + V (r )
σ

a.s.,

where E ∗ represents the expectation with respect to the bootstrap conditional distribution and V(r) is a
k × 1 Brownian motion with covariance matrix σ−2 I .

Thus the limiting distribution of a single bootstrap draw for the partial sums process mimics the result in
Proposition 3.4, except that the Brownian bridge UB (r ) is replaced with a Brownian motion V (r ). DefineY ∗ (r ) = r Y + V (r ). Using Proposition 3.7, we can obtain asymptotic representations for the different
procedures incorporating a bagging step in analogy with (3.7)-(3.11):

Proposition 3.8 In large samples, the distributions of the alternative parameter estimates including abagging step are as follows:(i) Using unrestricted MLE:
T 1/2 β̃ →d Σ−1 Y (1);15

(3.18)

(ii) Using the positive-part James-Stein estimator:

T

1/2

·
½
−1 ∗
β̃ →d E Σ Y (1) max 1 −
∗

¾¸
k −2
,0 ;
Y ∗ (1)0 Σ−2 Y ∗ (1)(3.19)(iii) Selecting the model using the AIC:

T 1/2 β̃ →d E ∗

·

X

¸
n
£
¤o
H (m ∗ )Y ∗ (1)1 m ∗ = arg max Y ∗ (1)0 H (m)Y ∗ (1) − 2n(m);mm∗(3.20)(iv) Selecting the model minimizing out-of-sample error:

T 1/2 β̃ →d E ∗hXn£
H (m ∗ )Y ∗ (1)1 m ∗ = arg max Y ∗ (1)0 H (m)Y ∗ (1)
m

m∗

¤oi
1
− Y ∗ (π)0 H (m)Y ∗ (π) + t r (H (m)Λ) log(π) ;π(3.21)(v) Using the split-sample method:T

1/2

β̃ →d E

∗

·

X

m∗

where

∗ Y

H (m )

∗

½
¸¾¸
·
1 ∗ 0
(1) − Y ∗ (π)∗∗
1 m = arg max
Y (π) H (m)Y (π) − 2n(m)
;
m
1−π
π

m ∗ denotes the summation over all the models in M , and E

P∗(3.22)

denotes the expectation with respect

to V (r ) in Y ∗ (r ) = r Y (1) + V (r ).

In Appendix B, we also provide more concrete expressions for the in-sample and split-sample procedures, in their standard form, with RB, and with bagging, in the special case where k = 1.
Bagging and our proposed RB procedure are closely related. RB uses simulation to integrate out the
Brownian bridge UB (r ). Bagging has an asymptotic representation that has a Brownian motion insteadof UB (r ), and then integrates out that Brownian motion. But RB has certain advantages. Bagging applies
only in the case of independent data, whereas our RB approach can be used in any setting where wehave consistent estimators of Λ and Σ. Also RB does not require resampling the data and reestimating
the model. This may make RB especially attractive when the data are dependent or in circumstanceswhere model estimation is computationally costly.
Breiman (1996) gave a heuristic argument for why bagging weakly reduces mean square error, but in fact

16

bagging can increase mean square error. Bühlmann and Yu (2002) showed this for the case of estimation
with AIC model selection. See also Andreas and Stuetzle (2000) and Friedman and Hall (2007). On theother hand RB is guaranteed to reduce the local asymptotic risk for any convex loss function.
Since RB involves integrating out UB (r ), it does not affect any procedure that does not depend on thisancillary noise. In particular, RB does not affect in-sample model selection with AIC. But bagging willaffect the limiting distribution of all of the procedures that we consider, except for the unrestricted MLE.

4 Extensions
In this section, we consider two extensions of the basic framework of our analysis, in the context of thelinear regression model.

4.1 Unmodeled Structural Change
A variant of our basic regression model specifies that y t = β0t x t + u t where T 1/2 β[Tr ] = W (r ), where W (r )may be either a stochastic or nonstochastic process. This allows various forms of structural breaks, andis similar to specifications used by Andrews (1993) and Elliott and Mueller (2014). For example, if βt =P
T −1/2 b + T −1/2 b̃1(t > [T s]), then W (r ) = b + b̃1(r > s). Or, if βt = T −1 ts=1 η s with Gaussian shocks, thenW (r ) is a Brownian motion. Proposition 4.1 gives the asymptotic distribution of the partial sum process
P ]T −1/2 σ−2 [Trt =1 x t y t in this variant of our basic model:

Proposition 4.1 As T → ∞, the partial sum process

T −1/2 σ−2[TrX]x t y t →d Z (r )

(4.1)

t =1
d Rr
where Z (r ) = Σ 0 W (s)d s +r ξ+UB (r ) , ξ ∼ N (0, Σ) and UB (r ) is an independent k-dimensional Brownian

bridge with covariance matrix Σ = σ−2 I .

Suppose that the researcher ignores the possibility of structural change, and simply uses the available
estimators for forecasting. The limiting distributions of the estimators will be as in Proposition 3.5, with17

Y (r ) replaced by Z (r ) everywhere. Alternatively, the researcher might be aware of the possibility of structural change, and might choose to select among models and estimate parameters using a rolling window.
Other approaches for dealing with the possibility of parameter instability might be considered, but weleave this topic for future research.
Structural instability may be part of the motivation for considering out-of-sample forecasting methods.
RB is no longer guaranteed to reduce risk in the above model with structural change. Nonetheless, theweb appendix also contains some simulations in which the parameters follow a random walk, but theeconometrician treats them as though they were fixed. Although there is no theoretical result that RB
has to reduce risk in this case, we find that risk is in fact lowered by RB as long as the parameter variation
is not too great.

4.2 Model Combination
It may also be appealing to combine forecasts made from multiple models, instead of selecting a single
model (Bates and Granger (1969) and Timmermann (2006)). Recalling that β̂(m) denotes the parameter
estimate from the model containing the variables indexed by m (with zeros in other locations), then we
could estimate the parameter vector as Σm w(m)β̂(m), where Σm denotes the sum over all the models
in M and the weights sum to 1. As examples of weighting schemes, we could set the weight for modelexp(AIC (m)/2)

m to w(m) = Σ ∗ exp(AIC (m ∗ )/2) (Buckland, Burnham, and Augustin, 1997) where AIC (m) denotes themAkaike Information Criterion in the model indexed by m, or we could weight models by out-of-sample
predictive performance setting:

i
`
(
β̂
(m))
t
1,t
−1
t =[T π]+1i.
hPw(m) = PT∗m ∗ expt =[T π]+1 `t (β̂1,t −1 (m ))

exp

hP
T

(4.2)

Alternatively, to do a combination version of the split-sample scheme, we could estimate the paramexp(AIC (π,m)/2)

eter vector as Σm w(m)β̂[T π]+1,T (m) where w(m) = P exp(AIC (π,m)/2) and AIC (π, m) denotes the AkaikemInformation Criterion for the model indexed by m computed only over the first fraction π of the sample.

Proposition 4.2 If the parameter vector is estimated by Σm w(m)β̂(m) then in large samples, the distribu-

18

tions of the alternative parameter estimates will be:X

w(m)H (m)Y (1),

(4.3)

¡
¢
w(m) ∝ exp [Y (1)0 H (m)Y (1) − 2n(m)]/2

(4.4)

¶
µ
1
0
0
w(m) ∝ exp Y (1) H (m)Y (1) − Y (1) H (m)Y (1) + t r (H (m)Λ) log(π) ,
π

(4.5)

m

where

or

for exponential AIC and out-of-sample prediction error weights, respectively. Meanwhile, if the parameter
vector is instead estimated by Σm w(m)β̂[T π]+1,T (m) with exponential AIC weights, then in large samples,the distribution of the estimator will be:X

Y (1) − Y (π)
,
1−π

(4.6)

¸ ¶
±
1
0
Y (π) H (m)Y (π) − 2n(m) 2 .
π

(4.7)

w(m)H (m)

m

where
µ·

w(m) ∝ exp

The standard bagging step can be added to any of these methods for forecast combination and the resulting limiting distributions are also given by Proposition 4.2, except with Y (·) replaced by Y ∗ (·) everywhere.
Or RB can be added, to integrate out UB (r ).
An alternative and more standard way to obtain combination weights for the out-of-sample forecasting
scheme would be to weight the forecasts by the inverse mean square error (Bates and Granger (1969)
and Timmermann (2006)). Under our local asymptotics, this will give each model equal weight in large
samples.

19

5 Numerical WorkIn this section we numerically explore the root mean squared errorq £¤E (T 1/2 β̃ − b)0 (T 1/2 β̃ − b) ,(5.1)

Given the expressions in Propositions 3.5 and 3.8, we can simulate the asymptotic risk of different methods in their standard form, with RB, and with bagging5 for different choices of the localization parameterb and the number of potential predictors k. None of the methods gives the lowest risk uniformly in b. Always using the big model is minmax, but due to the Stein phenomenon, it may be dominated by shrinkage estimators. In all cases, RB and bagging are implemented using 100 replications, the out-of-sample
and split-sample methods both set π = 0.5, and we set Σ = Λ = I k . The asymptotic risk is symmetric inb and is consequently shown only for non-negative b. The bagging results from Proposition 3.8 apply
only in the special case of the linear regression model, but RB applies in the general pseudo-likelihood
framework. Figure 1 plots the asymptotic risk of the standard MLE, JS, in-sample, out-of-sample andsplit-sample methods, for the case k = 1 against b. Results with RB and bagging are also included. The
web appendix includes risk calculations for the Laplace estimator (Magnus, 2002; Magnus, Powell, andPrüfer, 2010) and the LASSO estimator (Tibshirani, 1996).
Among the standard methods, selecting the model in-sample by AIC does better than the out-of-sample
scheme for most values of b, which in turn dominates the split-sample method. But RB reduces therisk of the out-of-sample and split-sample methods for all values of b, and makes them much morecompetitive. Bagging accomplishes much the same thing. The fact that bagging improves the out-ofsample and split-sample methods uniformly in b is just a numerical result, but it is also a theoreticalresult for RB. Neither bagging nor RB dominates the other in terms of risk. Bagging also helps with thein-sample method for some but not all values of b—this was also shown by Bühlmann and Yu (2002).
Among all the prediction methods represented in Figure 1, which one the researcher would ultimately
would want to use depends on b, which is in turn not consistently estimable. But the split-sample andout-of-sample methods do best for many values of the localization parameter, as long as the bagging or
RB step is included. This is especially true for small values of b.
5 The results with bagging are based on Proposition 3.8, which applies only in the case of the linear regression model.

20

We next consider the case where the number of potential predictors k is larger, but only one parameter
actually takes on a nonzero value. (Of course, the researcher does not know this.) Without loss of generality, we let the nonzero element of b be the first element and so specify that b = (b 1 , 0, . . . , 0)0 . To keep themodel space manageable, the set of possible models M consists of models with the first l predictors, forl ∈ {0, 1, ...k}. Figures 2 and 3 plot the risk for k = 3 and k = 20 against b 1 for in-sample, out-of-sample andsplit-sample methods in the standard form, with RB, and with bagging. The positive-part James-Steinestimator is also included. The split-sample method with either RB or bagging compares very favorablywith the other alternatives. Other numerical results with multiple predictors are contained in the web
appendix, and give similar conclusions.
Figure 4 plots the risk for k = 1 against b for the in-sample, out-of-sample and split-sample forecast
combination methods, in their standard form, with RB and with bagging. These are based on simulating
the distributions in Proposition 4.2. The combination forecasts are generally better than forecasts basedon selecting an individual model. Nonetheless, with combined forecasts as with individual forecasts, inthe absence of a randomization step, using in-sample AIC weights does best for most values of b. But
adding in RB/bagging reduces the risk of the combination forecasts with out-of-sample or split-sample
weights uniformly in b, and makes these methods substantially more competitive.
In the web appendix, we also carry out small-sample Monte Carlo simulations. These give similar conclusions to the local asymptotic calculations reported here, although estimation of the variance-covariancematrices does increase MSPE when k is large.

6 Application
We finally consider a small empirical illustration to the classic problem of forecasting stock returns. We
follow the setup of Goyal and Welch (2008) in forecasting annual excess stock returns using data from1926 to 2014 with 13 possible predictors: book-to-market ratio, Treasury bill yields, long-term yields,
net equity expansion, inflation, percent equity issuing, long term returns, stock variance, default yieldspread, default return spread, dividend-price ratio, dividend yield and earnings-price ratio. These are
all the predictors considered by Goyal and Welch (2008) for which data are available over the full sample
period, excepting some that are perfectly multicollinear with the included predictors. We then consider
21

one-year-ahead forecasting of stock returns using all possible subsets of these 13 predictors including
the empty set, for a total of 8,192 models. Definitions of the predictors and sources are in Goyal andWelch (2008). Each model includes an intercept.
AIC chooses a single predictor: the book-to-market ratio. The same is true selecting the model with the
split-sample approach. Selecting the model by out-of-sample performance, percent equity issuing andthe dividend price ratio are the two chosen predictors. Table 1 shows the probability of the chosen model
having no predictors whatsoever (other than the intercept), using AIC, out-of-sample and split-sample
schemes with bagging and out-of-sample and split-sample schemes with RB. Having no predictors otherthan the intercept amounts to forecasting excess stock returns using the unconditional mean. Table 1
also shows the expected number of predictors (other than the intercept) for all five simulation schemes.
AIC with bagging, AIC with RB and the split-sample approach with RB all put most weight on models
with 0 or1 predictors, and thereby implying heavy shrinkage. The out-of-sample approach with either
bagging or RB chooses larger models. Table 1 also reports the fraction of draws for which each of the
13 possible predictors is included in the model, for each of the five simulation schemes. Finally, Table 1
gives the forecast for excess stock returns in 2015 using each of the methods. These are very close to theunconditional mean, indicating a considerable degree of shrinkage.

7 Conclusion
When forecasting using k potential predictors, each of which has a coefficient that is local to zero, there
are several competing methods, none of which is most accurate uniformly in the localization parameter,
which is in turn not consistently estimable. Optimizing the in-sample fit, as measured by the Akaike
information criterion, generally does better than out-of-sample or split-sample methods. However, theout-of-sample and split-sample methods can be improved substantially by removing the impact of anancillary noise term that appears in their limit representations, either through Rao-Blackwellization or
bagging. Rao-Blackwellization uniformly lowers asymptotic risk of the out-of-sample and split-sample
methods. For important ranges of the local parameters, these modified procedures are very competitive
with in-sample methods.

22

A Appendix: Proofs of Propositions
Proof of Proposition 3.4: We haveT −1/2[Tr
X] ∂`t (0)
t =1

∂β

= T −1/2[TrX] ∂`t (β0,T )
t =1

∂β

− T −1/2[TrX] ∂2 `t (β0,T )β0,T + o p (1)
∂β∂β0
t =1

→d r Σb + B (r ),where B (r ) denotes a Brownian motion with covariance matrix Λ. Let Y (r ) = r Σb+B (r ), and let Y = Y (1)which is N (Σb, Λ). Define UB (r ) = B (r ) − r B (1). Then UB is a Brownian bridge by standard arguments,
and by calculating the covariance between UB (r ) and B (1) it can be verified that UB is independent of
B (1). We can therefore write

Y (r ) = r Σb + B (r ) = r Σb + r B (1) +UB (r )= r Y +UB (r ).

and Y = Σb + B (1) is uncorrelated with UB , and hence independent of UB . ■Proof of Proposition 3.5: Let β̂ denote the unrestricted estimator and let β̂(m) denote the restricted
estimators as defined in Section 2. Equations (3.7) and (3.8) immediately follow because T 1/2 β̂ →d Y .
The AIC objective function (to be maximized) is:

2

T
X
t =1

l t (β̂(m)) − 2n(m) = 2

T
X
t =1

l t (0) + 2β̂(m)0

T ∂l (0)
X
t
t =1

∂β

+ β̂(m)0

T ∂2 l (0)Xtt =1 ∂β∂β

0β̂(m) − 2n(m) + o p (1),which is asymptotically the same, up to the same affine transformation across all models, asY (1)0 H (m)Y (1) − 2n(m),noting that H (m)ΣH (m) = H (m).

23

The OOS objective function (to be maximized) is:·¸2
1
0 ∂ l t (0)0 ∂l t (0)l t (0) + β̂1,t −1 (m)l t (β̂1,t −1 (m)) =+ β̂1,t −1 (m)β̂1,t −1 (m) + o p (1).
∂β
2
∂β∂β0
t =[T π]+1t =[T π]+1
T
X

T
X

Let S(m) be the diagonal matrix with ones in the diagonal elements indexed by m and zeros elsewhere.
Now T 1/2 β̂1,[Tr ] (m) →d H (m) Y r(r ) andTXt =[T π]+1

β̂1,t −1 (m)00tX
−1 ∂l (βTX
∂l t (β0,T )∂l t (0)1s 0,T )
=
H (m)
∂β
∂β
t
−
1
∂β
t =[T π]+1 s=1+ b 0 ΣH (m)T −1/2TX

∂l t (β0,T )t =[T π]+1

∂β

+

T
X
t =[T π]+1

β̂1,t −1 (m)01
Σb + o p (1)
T

0

Z 1
B (r )Y (r ) 0→d
H (m)Σbd r
H (m)d B (r ) + b 0 ΣH (m)(B (1) − B (π)) +
r
r
π
π
Z 1
Z 1
Y (r )0
Y (r ) 0=
H (m)d B (r ) +H (m)Σbd r
r
r
π
π
Z 1
Y (r )0=
H (m)d Y (r ).
rπZ 1Consequently, the OOS objective function is asymptotically the same, up to the same affine transformation across all models, as
Z 1
Y (r )0
Y (r )0
Y (r )
2
H (m)d Y (r ) −
H (m)
d r.
r
r
r
π
π
Z 1

Using an argument similar to Hansen and Timmermann (2015), define1
F (Y , r ) = Y (r )0 H (m)Y (r ) − t r (H (m)Λ) log(r ).
r
By Ito’s lemma
1
1
1
2d F (Y , r ) = {− 2 Y (r )0 H (m)Y (r ) − t r (H (m)Λ) + 2b 0 ΣH (m)Y (r ) + t r (H (m)Λ)}d r + Y (r )0 H (m)d B (r )
r
r
r
r
1
2
0
0
0
= {− 2 Y (r ) H (m)Y (r ) + 2b ΣH (m)Y (r )}d r + Y (r ) H (m)d B (r )rr
12= {− 2 Y (r )0 H (m)Y (r ) + 2b 0 ΣH (m)Y (r )}d r + Y (r )0 H (m)d Y (r ) − 2Y (r )0 H (m)Σbd r
rr
2
1= Y (r )0 H (m)d Y (r ) − 2 Y (r )0 H (m)Y (r )d r.
r
r

24

So
Z 1
Y (r )0
Y (r )0
Y (r )
H (m)d Y (r ) −H (m)drF (Y , 1) − F (Y , π) =d F (Y , r ) = 2
r
r
r
π
π
π
1
= Y (1)0 H (m)Y (1) − Y (π)0 H (m)Y (π) + t r (H (m)Λ) log(π).
π
Z 1

Z 1

Consequently, the OOS objective function is asymptotically the same, up to the same affine transformation across all models, asY (1)0 H (m)Y (1) − π1 Y (π)0 H (m)Y (π) + t r (H (m)Λ) log(π).
The AIC estimated over the first fraction π of the sample is:

2

[T
π]
X

l t (β̂1,[T π] (m)) − 2n(m) = 2

t =1

[T
π]
X

l t (0) + 2β̂1,[T π] (m)0

t =1

[T
π] ∂l (y , 0)
X
t t
t =1

+ β̂1,[T π] (m)0∂β[Tπ] ∂2 l (y , 0)Xt tβ̂1,[T π] (m) − 2n(m) + o p (1)
0
∂β∂β
t =1

which is asymptotically the same, up to the same affine transformation across all models, as

10
π Y (π) H (m)Y (π) − 2n(m)The limiting distributions in Proposition 3.5 all follow from these results and the facts that T 1/2 β̂ →dΣ−1 Y (1) and T 1/2 β̂(m) →d H (m)Y (1).

■Proof of Proposition 3.6: We know that T 1/2 β̂(m) →d H (m)Y (1) = H (m)Ỹi (1) and Ĥ (m) →p H (m).
U i (r )Hence T 1/2 β̂(m) + Ĥ (m) Br→d H (m) Ỹir(r ) . The result follows immediately.

■

Proof of Proposition 3.7: Let {x t∗ , y t∗ } denote a bootstrap sample and let u t∗ = y t∗ −β00,T x t∗ , t = 1, . . . , T .β̂ isthe OLS estimate using the original data. Using the assumption that the x t are orthonormal, we can takePβ̂ = T −1 Tt=1 x t y t . Then, in the bootstrap conditional probability space (conditioning on the realization25

of {x t , y t }),TX] ∗ ∗ r p
X] ∗ ∗
X
1 [Tr
1 [Tr
r
x t y t − 2 T β̂ = p
xt yt − pxt ytpσσ2 T t =1
σ2 T t =1
σ2 T t =1
TX] ∗ ¡ ∗0
X
¢
1 [Tr
r
= p
x t x t β0,T + u t∗ − p
x t (x t0 β0,T + u t )σ2 T t =1
σ2 T t =1
"
Ã
!
#!
Ã
T
Tr
T
X] ∗ ∗0 1 X
1 X
1 X
1 [Tr
0
∗ ∗
xt xt −x s x s β0,T + p
xt ut −
xs us
=
p
T s=1
T s=1
σ2 T t =1
σ2 T t =1
Ã
Ã
"
!
#!
Tr
T
T
X] ∗ ∗0 1 X
1 X
1 [Tr
1 X
∗ ∗
0
xt ut −
= 2x x −xs xs b + p
xs us .
σ T t =1 t t
T s=1
T s=1
σ2 T t =1

The first term is o p∗ (1) by a Uniform Law of Large Numbers for partial sums (e.g., Gaenssler and Ziegler(1994)), and the second converges almost surely to V (r ), from Theorem 2.2 of Park (2002). Therefore
X] ∗ ∗
1 [Trr p
x t y t = 2 T β̂ + V (r ) + o p (1) a.s.
p
σ
σ2 T t =1

Given that g is uniformly integrable,
" Ã

E

∗

X] ∗ ∗
1 [Tr
g
xt yt
p
σ2 T t =1!#

h ³ r p
´i
→ E g 2 T β̂ + V (r )
σ

a.s.,

as required.

■

The proofs of Propositions 3.8 and 4.2 involve exactly the same calculations as in Proposition 3.5 and arehence omitted.

Proof of Proposition 4.1: We haveT −1/2 σ−2[TrX]x t y t = T −1/2 σ−2

t =1

[Tr
X]

x t x t0 βt + T −1/2 σ−2

t =1= T −3/2 σ−2

[Tr
X]
t =1

→d ση Σ

Z r
0

[Tr
X]

xt ut

t =1

x t x t0 Σts=1 η s + T −1/2 σ−2W (s)d s + B (r ) = ση Σwhere B (r ) is a Brownian motion with covariance matrix Λ.
26

Z r
0

[Tr
X]

xt ut

t =1

W (s)d s + r ξ +UB (r ).

■B Appendix: Shrinkage Representations in the case k = 1
In the case k = 1, and with Σ = Λ, some of the expressions in Propositions 3.5 and 3.8 can be simplified.
For the AIC estimator in its standard form we have:³
p ´
T 1/2 β̃ →d Σ−1 Y 1 |Y | > 2Σ .

(B1)

For the AIC estimator, with bagging, we have:
o
n
p
p
p
p
T 1/2 β̃ →d Σ−1 Y − Y Φ( 2 − κY ) + κφ( 2 − κY ) + Y Φ(− 2 − κY ) − κφ(− 2 − κY ) ,where κ = Σ−1/2 , shown in proposition 2.2 of Bühlmann and Yu (2002).6 Comparing this to (B1), bagging
smooths the hard thresholding in the original AIC estimator.
For the split-sample estimator, we have:sÃ

T 1/2 β̃ →d Σ−1 z 1 1 |z 2 | >

!
2Σ
,
π

1
B (π)
and z 2 = Y + UBπ(π) . By direct calculations, z 1 is N (Σb, 1−πΣ), z 2 is N (Σb, π1 Σ) and z 1
where z 1 = Y − U1−π

and z 2 are mutually independent. With RB, for the split-sample estimator in the i th simulated sample,
we have:
s
s !
! Ã
πΣ
(1
−
π)Σ
2Σ
T 1/2 β̃i →d Σ−1 Y −z(i ) 1 |Y +z(i )| >
1−π
π
π
Ã
!
r
2
= (Σ−1 Y − γz(i ))1 |γY + z(i )| >
,
1−π
Ã

s

where z(i ) is N (0, 1), and is independent of Y and γ =qπ(1−π)Σ . Thus for the overall split-sample estima-

6 Indeed, given the orthonormal setting, even if k > 1, if we sort the coefficient estimates by their absolute magnitude andapply AIC sequentially to these models, dropping variables one at a time as long as called for by the information criterion, thenthe above two expressions will apply to each element of β̃ − β (Bühlmann and Yu, 2002; Stock and Watson, 2012). But the use of
the AIC that we are considering in this paper is to select among all 2k possible models and so no such simplification is available
in this case.

27

tor with RB, we have:
(

T

1/2

−1

β̃ →d ΣÃrY −Y Φ

!
Ãr
!
Ã r
!
Ãr
!)
2
2
2
2
− γY − γφ− γY + Y Φ −− γY + γφ
− γY
.
1−π
1−π
1−π
1−π

Meanwhile, for bagging the split-sample estimator in the i th bagging sample, we have:sÃ

T 1/2 β̃i →d z 1 (i )1 |z 2 (i )| >

2Σ
π

!

Σ
i (π)
and z 2 (i ) = Y + Viπ(π) . By direct calculations, z 1 (i )|Y is N (Y , (1−π)), z 2 (i )|Y iswhere z 1 (i ) = Y + Vi (1)−V
1−π

N (Y , Σπ ) and the two are independent, conditional on Y . Thus the split-sample with bagging estimator is:

T

1/2

−1

β̃ →d Σ

r
r
½
µ
¶
µ
¶¾
p
p
π
π
Y −Y Φ 2−
Y +Y Φ − 2−Y .
ΣΣFor the out-of-sample estimator, we have:¶µq
T 1/2 β̃ →d Σ−1 Y 1 |z 3 | < max(Y 2 π + π log(π)Σ, 0) ,where z 3 = πY + UB (π) is N (πΣb, πΣ). With RB, for the out-of-sample estimator, in the i th simulated
sample, we have:Ã

T 1/2 β̃i →d Σ−1 Y 1 |z(i )| <s!
¶log(π)
Y2+max,0 −b .
πΣ2πΣµThus for the overall out-of-sample estimator with RB:
!
Ã s
!)
¶
¶
µ 2
Y2
log(π)Ylog(π)
max+, 0 − b − Φ − max+,0 −b .
πΣ2
πΣ
πΣ2
πΣ

( Ãs

T 1/2 β̃ →d Σ−1 YΦµLastly, for bagging the out-of-sample estimator in the i th bagging sample, we have:©
ª
T 1/2 β̃i →d Σ−1 (Y + Vi (1))1 (Y + Vi (1))2 − (πY + Vi (π))2 + log(π)Σ > 0 .

(We do not have a closed form expression for the overall out-of-sample estimator with bagging.)
Thus in all cases we can write the limit of T 1/2 β̃ as Y g (Y ).

28

References
A KAIKE , H. (1974): “A New Look at the Statistical Model Identification,” IEEE Transactions on AutomaticControl, 19, 716–723.
A NDREAS , B., AND W. S TUETZLE (2000): “Bagging does not always Decrease Mean Squared Error,” mimeo.
A NDREWS , D. W. K. (1993): “Tests for Parameter Instability and Structural Change with Unknown Change
Point,” Econometrica, 61(4), 821–856.
A SHLEY, R., C. W. G RANGER , AND R. S CHMALENSEE (1980): “Advertising and Aggregate Consumption: An
Analysis of Causality,” Econometrica, 48, 1149–1167.
B ATES , J. M., AND C. W. G RANGER (1969): “The combination of forecasts,” Operations Research Quarterly,
20, 451–468.
B REIMAN , L. (1996): “Bagging Predictors,” Machine Learning, 36, 105–139.
B UCKLAND, S. T., K. P. B URNHAM , AND N. H. AUGUSTIN (1997): “Model Selection: An Integral Part ofInference,” Biometrics, 53, 603–618.
B ÜHLMANN , P., AND B. Y U (2002): “Analyzing Bagging,” Annals of Statistics, 30, 927–961.
C HRISTOFFERSEN , P., AND F. X. D IEBOLD (1997): “Optimal Prediction Under Asymmetrical Loss,” Econometric Theory, 13, 806–817.
C LAESKENS , G., AND N. L. H JORT (2008): Model Selection and Model Averaging. Cambridge University
Press, Cambridge.
C LARK , T. E. (2004): “Can Out-of-Sample Forecast Comparisons help Prevent Overfitting?,” Journal of
Forecasting, 23, 115–139.
D AVIDSON , J. (1994): Stochastic Limit Theory. Oxford University Press.
D IEBOLD, F. X., AND R. S. M ARIANO (1995): “Comparing Predictive Accuracy,” Journal of Business andEconomic Statistics, 13, 253–263.
E FRON , B. (2014): “Estimation and Accuracy After Model Selection (with discussion),” Journal of theAmerican Statistical Association, 109, 991–1007.
29

E LLIOTT, G., AND U. K. M UELLER (2014): “Pre and Post Break Parameter Inference,” Journal of Econometrics, 180(2), 141–157.
F RIEDMAN , J. H., AND P. H ALL (2007): “On Bagging and Nonlinear Estimation,” Journal of Statistical
Planning and Inference, 137, 669–683.
G AENSSLER , P., AND K. Z IEGLER (1994): “A Uniform Law of Large Numbers for Set-Index Processes
with Applications to Empirical and Partial-Sum Processes,” in Probability in Banach Spaces, 9, ed. by
J. Hoffmann-Jorgensen, J. Kuelbs, and M. B. Marcus. Springer.
G IACOMINI , R., AND H. W HITE (2006): “Tests of Conditional Predictive Ability,” Econometrica, 74, 1545–
1578.
G OYAL , A., AND I. W ELCH (2008): “A Comprehensive Look at The Empirical Performance of Equity Premium Prediction,” Review of Financial Studies, 21, 1455–1508.
H ANSEN , P. R. (2009): “In-Sample Fit and Out-of-Sample Fit: Their Joint Distribution and its Implicationsfor Model Selection,” mimeo.
H ANSEN , P. R., AND A. T IMMERMANN (2015): “Equivalence Between Out-of-Sample Forecast Comparisons and Wald Statistics,” Econometrica, 83(6), 2485–2505.
I NOUE , A., AND L. K ILIAN (2004): “In-Sample or Out-of-Sample Tests of Predictability: Which One Should
We Use?,” Econometric Reviews, 23, 371–402.
(2006): “On the Selection of Forecasting Models,” Journal of Econometrics, 130, 273–306.
L EEB , H., AND B. M. P ÖTSCHER (2005): “Model selection and inference: Facts and Fiction,” Econometric
Theory, 21, 21–59.
M AGNUS , J. R. (2002): “Estimation of the Mean of a Univariate Normal Distribution with Known Variance,” Econometrics Journal, 5, 225–236.
M AGNUS , J. R., O. P OWELL , AND P. P RÜFER (2010): “A Comparison of Two Model Averaging Techniqueswith an Application to Growth Empirics,” Journal of Econometrics, 154, 139–153.
M ALLOWS , C. L. (1973): “Some Comments on Cp,” Technometrics, 15, 661–675.
30

M EESE , R. A., AND K. R OGOFF (1983): “Empirical Exchange Rate Models of the Seventies: Do They Fit
Out of Sample?,” Journal of International Economics, 14, 3–24.
PARK , J. (2002): “An Invariance Principle for Sieve Bootstrap in Time Series,” Econometric Theory, 18,
469–490.
S TOCK , J. H., AND M. W. WATSON (2012): “Generalized Shrinkage Methods for Forecasting Using Many
Predictors,” Journal of Business and Economic Statistics, 30, 481–493.
T IBSHIRANI , R. (1996): “Regression Shrinkage and Selection via the Lasso,” Journal of the Royal StatisticalSociety, Series B, 58, 267–288.
T IMMERMANN , A. (2006): “Forecast Combination,” in Handbook of Economic Forecasting, ed. by C. W.
Granger, G. Elliott, and A. Timmermann, Amsterdam. North Holland.
VAN D ER VAART, A. W. (1998): Asymptotic Statistics. Cambridge University Press, Cambridge.
W EST, K. D. (2006): “Forecast Evaluation,” in Handbook of Economic Forecasting, ed. by C. W. Granger,
G. Elliott, and A. Timmermann, Amsterdam. North Holland.
W ILSON , E. (1934): “The Periodogram of American Business Activity,” 34, 375–417.

31

Table 1: Application to Goyal-Welch Data

AICB OOSBProbability of no predictors
0.72
0.00
Expected Number of predictors
0.28
4.86
Marginal probability of inclusion of each predictor:
Book-to-market ratio
0.20
0.40
Treasury bill yields
0.08
0.36
Long-term yields
0.00
0.20
Net equity expansion
0.00
0.34
Inflation
0.00
0.32
Percent equity issuing
0.00
0.78
Long term returns
0.00
0.54
Stock variance
0.00
0.32
Default yield spread
0.00
0.40
Default return spread
0.00
0.26
Dividend-price ratio
0.00
0.20
Dividend yield
0.00
0.30
Earnings-price ratio
0.00
0.44
Forecast for Excess Returns in 2015 0.05
0.05
Memo: Unconditional Mean, 1926-2014: 0.06

SSB
0.40
0.60

OOSRB
0.00
5.00

SSRB
0.34
0.66

0.30
0.14
0.00
0.14
0.00
0.02
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.05

0.46
0.68
0.48
0.08
0.22
0.86
0.80
0.42
0.52
0.00
0.12
0.20
0.16
0.08

0.58
0.02
0.00
0.06
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.04

Notes: For the empirical application described in Section 7, this table reports the fraction of draws witheach simulation scheme that includes no predictors (other than the intercept), the average number ofpredictors (other than the intercept) chosen in each scheme, and the fraction of draws that includes eachone of the predictors. The objective is forecasting year-ahead excess stock returns using annual data from
1926-2014. The forecast for excess stock returns in 2015 and the unconditional mean are also included.
Model selection is over all 8,192 permutations of the 13 predictors. Definitions of the predictors and data
sources are in Goyal and Welch (2008).

32

Figure 1: Local Asymptotic Risk (k = 1)
No resampling

Bagging

2

2

1.5

1.5

1

1

0.5

0.5
0

2

4

6

0

b

2

4

6

b

Rao-Blackwellization
2

1.5

AIC
OOS
Split Sample
MLE

1

0.5
0

2

4

6

b
Notes: These are the simulated local asymptotic risk values, equation (5.1), for different procedures, plotted against b.
Note that MLE is the same without any resampling, with bagging or with RB. AIC is the same without any resampling orwith RB.

33

Figure 2: Local Asymptotic Risk (k = 3)
No resampling

Bagging

2

2

1.5

1.5

1

1

0.5

0.5
0

2

4

6

0

b

2

4

6

b

1

1

Rao-Blackwellization
2
James Stein
AIC
OOS
Split Sample
MLE

1.5

1

0.5
0

2

4

6

b1
Notes: These are the simulated local asymptotic risk values, equation (5.1), for different procedures, plotted against
b 1 , where b = (b 1 , 0, . . . , 0)0 . Note that MLE is the same without any resampling, with bagging or with RB. AIC is the same
without any resampling or with RB.

34

Figure 3: Local Asymptotic Risk (k = 20)
No resampling

Bagging

4.5

4.5

4

4

3.5

3.5

3

3

2.5

2.5

2

2

1.5

1.5

1

1
0

2

4

6

0

b

2

4

6

b

1

1

Rao-Blackwellization
4.5
4
3.5

James Stein
AIC
OOS
Split Sample
MLE

3
2.5
2
1.5
1
0

2

4

6

b1
Notes: These are the simulated local asymptotic risk values, equation (5.1), for different procedures, plotted against
b 1 , where b = (b 1 , 0, . . . , 0)0 . Note that MLE is the same without any resampling, with bagging or with RB. AIC is the same
without any resampling or with RB.

35

Figure 4: Local Asymptotic Risk: Combination Forecasts (k = 1)
No resampling

Bagging

1.6

1.6

1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4
0

2

4

6

0

b

2

4

6

b

Rao-Blackwellization
1.6
1.4
Exponential AIC Weights
Exponential OOS Weights
Exponential Split Sample Weights
MLE (No Combination)

1.2
1
0.8
0.6
0.4
0

2

4

6

b
Notes: These are the simulated local asymptotic risk values, equation (5.1), for different procedures, plotted against b.
Note that MLE is the same without any resampling, with bagging or with RB. Exponential AIC forecast combination is thesame without any resampling or with Rao-Blackwellization.

36