title: Minimizing post-shock forecasting error through aggregation of outside information
authors: Jilei Lin Daniel J. Eck
text: 
  
   The polar-generalized normal distribution
      
   Masoud  Faridi,
  Majid Jafari Khaledi^a
  
  masoudfaridi@modares.ac.ir
  
  ^a Department of Statistics, Tarbiat Modares University, P.O. Box 14115-134, Tehran, Iran
  
  jafari-m@modares.ac.ir
  
      
  ============================================================================================================================================================================================
  
  
      We develop a forecasting methodology for providing credible forecasts for time series that have recently undergone a shock. We achieve this by borrowing knowledge from other time series that have undergone similar shocks for which post-shock outcomes are observed. Three shock effect estimators are motivated with the aim of minimizing average forecast risk. We propose risk-reduction propositions that provide conditions that establish when our methodology works. Bootstrap and leave-one-out cross validation procedures are provided to prospectively assess the performance of our methodology. Several simulated data examples, and a real data example of forecasting Conoco Phillips stock price are provided for verification and illustration.
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  We provide forecasting adjustment techniques with the goal of lowering overall forecast error when the time series under study has undergone a structural shock. We focus on the specific setting in which a structural shock has occurred and one desires a prediction for the post-shock response at the next time point. Standard forecasting methods may not yield accurate predictions in the presence of such structural shocks <cit.>. This is a general problem that has many real life applications. For example, one may be interested in forecasting the stock price of a company tomorrow after hearing terrible or great news about the company  after hours trading. Companies may be interested in forecasting the demand of their products to adjust production after they were involved in a brand crisis, but they only have recent sales data for which the company is operating well. All is not lost in this setting, one may be able to supplement the present forecast with past data borrowed from other time series which contain similar structural shocks. The core idea of our methodology is to sensibly aggregate similar past realized shock effects which arose from other time series, and then incorporate the aggregated shock effect estimator into the present forecast. Our method of combining shock effects embraces ideas from conditional forecasting <cit.>, time series pooling using cross-sectional panel data <cit.>, forecasting with judgement and models <cit.>, synthetic control methodology <cit.>, and expectation shocks <cit.>. 
  
  We study the post-shock forecasting problem in the context of additive shock effects in linear autoregressive models. In this post-shock forecasting setting, the researcher has a time series of interest which is known to have recently undergone a structural shock, and the post-shock response is not observed. In this setting, the additive shock effect is a random effect that is parameterized in the autoregressive model. The shock effect is then estimated using ordinary least squares (OLS). The researcher must move beyond the modeling paradigm that they were previously working under to accommodate this new shock effect <cit.>. One method for estimating the shock effect is to produce a conditional point forecast where a sequence of non-zero future structural shocks are conditioned upon and estimated <cit.>. Such conditional point forecasts are appropriate when the shock sequence considered is within the range of historical experience <cit.>. On the other hand, our methodology allows for the inclusion of outside data sources and covariates into this conditional forecasting context provided that the shock effects from outside data sources are all thought to arise from a  data generating process similar to that of the shock under study. For further differences of assumptions on shocks, our methodology allows for unprecedented shocks and no observation of past shocks. 
  
  In our methodological framework, the researcher creates a synthetic panel of time series which have undergone similar structural shocks in the past. Construction of the donor pool that forms this synthetic panel is similar to that in synthetic control methodology (SCM) <cit.>. As in SCM, care is needed when forming the donor pool of time series. However, there are key differences between our framework and SCM. We assume that the time series forming the donor pool are independent from the time series under study before the timing of the shock. Moreover, the shocks from candidate time series in the donor pool, together with the shock in the time series of interest, are assumed to be  from a common family of distributions with existing first and second moments. 
  
  We estimate the shock effects that are present in the time series forming the donor pool for which post-shock responses are observed. We then aggregate these estimated shock effects and use this aggregated estimate as an estimator for the shock effect in the time series of interest. This estimator is then added to a forecast for the yet to be realized post-shock response corresponding to the time series of interest. Shock effects in our post-shock forecasting framework is similar to “expectation shocks” which are studied in <cit.>. The context in <cit.> allows for consistent estimation of expectation shocks under a vector autoregressive model, possibly involving an instrumental variable approach as in <cit.>. In our context, the yet to be observed shock effect of interest is a random effect, and we can only partially estimate features of the random effect distribution using the time series forming the donor pool. 
  
  
  
  
  
  
  
  
  
  
  In this article, we will assume a simple autoregressive data generating process similar to that in <cit.> with a general random-effect structure. Therefore, our methodology is similar to the “K latent pooling” framework of <cit.>. However, our model formulation is more general than <cit.>. In our model, the donor pool can consist of dependent time series but time series within the donor pool should be independent of the time series of interest. However, mutual independence among time series in the donor pool can aid prospective evaluation of the reliability of our method. We consider three aggregation techniques: simple averaging, inverse-variance weighted averaging, and similarity weighting. The latter technique is similar to the weighting in synthetic control methodology <cit.>. Our autoregressive model will consider present day covariates to better motivate similarity weighting. The considered adjustment strategies all target the mean of the shock effect distribution. Such an estimation strategy can reduce mean squared error (MSE) when variation in the shock effect distribution is small relative to the mean. We provide risk-reduction propositions that detail the conditions when the adjusted forecasts will work better than the original forecast. The involved parameters in the risk-reduction propositions can be estimated by a residual bootstrap procedure that we develop. We also motivate a simple leave-one-out cross validation procedure which can prospectively assess the performance of our shock effect adjustment estimators. This prospective assessment does not require the observation of the post-shock response. Our Monte Carlo simulation results show that the risk-reduction propositions are nearly perfectly correct when the model for the shock effects is identified well with appropriate covariates under a fixed design. We demonstrate the utility of our methodology in a real data analysis in which we forecast the stock price of Conoco Phillips shares that experienced a large structural shock on March 9th, 2020. We will show that our proposed adjustment estimators yield much better results than no adjustment in this setting. We also use this example to demonstrate settings in which the shock effect may be decomposed into separate estimable parts. We now motivate our framework for post-shock forecasting.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § SETTING
  
  
  
  We will suppose that a researcher has time series data (y_i,t,_i,t), for t = 1, …,  T_i and i = 1, …, n+1, where y_i,t is a scalar response and _i,t is a vector of covariates that are revealed to the analyst prior to the observation of y_1,t.  Suppose that the analyst is interested in forecasting y_1,t, the first time series in the collection. We will suppose that specific interest is in forecasting the response after the occurrence of a structural shock. To gauge the performance of forecasts, we consider forecast risk in the form of MSE,
  
      R_T = 1/T∑_t=1^Tŷ_1,t - y_1,t^2,
  
  and root mean squared error (RMSE), given by √(R_T), in our analyses. In this article, we focus on post-shock prediction where forecasts methods only differ at the next future time point. Thus the MSE reduces to the magnitude ŷ_1,t - y_1,t^2.
  
  
  
  Our post-shock forecasting methodology will consist of selecting covariates _i,t, constructing a suitable donor pool of candidate time series that have undergone similar structural shocks to the time series under study, and specifying a model for the time series (y_i,t,_i,t), for t = 1, …,  T_i and i = 1, …, n+1. In this article, we consider a dynamic panel data model with autoregressive structure similar to that in <cit.>. Our dynamic panel model includes an additional shock effect whose presence or absence is given by the binary variable D_i,t, and we will assume that the donor pool time series are independent of the time series under study. The details of this model are in the next section.
  
  
  
  Figure <ref> provides a simple intuition of the practical usefulness of our proposed methodology. This figure depicts a time series that experienced a shock at time point T_1^*+1 = 126. It is supposed that the researcher does not have any information beyond T_1^*+1, but does have observations from a donor pool of forty time series that have previously undergone a similar shock for which post-shock responses are recorded. Similarity in this context means that the shock effects are random variables that from a common distribution. In this example, the mean of the estimated shock effects is taken as a shock effect estimator for the time series under study. Forecasts are then made by adding this shock effect estimator to the estimated response values obtained from the estimation procedure that ignores the shock. It is apparent from Figure <ref> that adjusting forecasts in this manner 1) leads to a reduction in forecasting risk; 2) does not fully recover the true shock effect. We evaluate the performance of this post-shock forecasting methodology throughout this article; we outline situations for when it is expected to work and when it is not.
  
  
  
  
  
  
  
  
  
  
   §.§ Model Setup
  
  
  
  
  In this section, we will describe the assumed dynamic panel models for which 
  post-shock aggregated estimators are provided. The basic structures of these models 
  are the same for all time-series in the analysis, the differences between them lie in the setup of the shock effect distribution.
  
  Let I(·) be an indicator function, T_i be the time length of the time series i for i = 1, …, n+1, and T_i^* be the time point just before the one when the shock is known to occur, with T_i^* < T_i.  For t= 1, …, T_i and i = 1, …, n+1, the model M_1 is defined as
  
      M_1  y_i,t =η_i +α_i D_i,t + ϕ_i y_i, t-1 + θ_i'x_i,t + ε_i,t
  
   where D_i,t = I(t = T_i^* + 1) 
  and _i,t∈^p with p ≥ 1.  We assume that the 
  x_i,t's are fixed. Let |x| denote the absolute value of x for x∈. For i = 1, …, n+1 and t=1, …, T_i, the random effects structure for M_1 is:
  
      η_i    F_η with  E_F_η(η_i) = 0, Var_F_η(η_i)  = σ^2_η
      ϕ_i    F_ϕ where  |F_ϕ| < 1, 
      θ_i    F_θ with  E_F_θ(θ_i) = μ_θ, Var_F_θ(θ_i)  = Σ^2_θ
      α_i    F_α with  E_F_α(α_i) = μ_α, Var_F_α(α_i)  = σ^2_α
      ε_i,t   F_ε with  E_F_ε(ε_i,t) = 0, Var_F_ε(ε_i,t)  = σ^2 σ > 0 ,  
      η_i    α_i ϕ_i θ_i ε_i,t.
  
  Notice that M_1 assumes that α_i are iid with α_i=μ_α 
  for i = 1, …, n+1. We also consider a model where the shock effects are linear functions of covariates with an additional additive mean-zero error. For i = 1, …, n+1, the random effects structure for this model (model M_2) is:
  
      M_2 [ y_i,t =η_i +α_i D_i,t + ϕ_i y_i, t-1 + θ_i'x_i,t + ε_i,t;                         α_i = μ_α+δ_i'x_i, T_i^*+1+ ε_i, ]
  
   where the added random effects are
  
      ε_i   F_ε with E_F_ε(ε_i)=0, Var_F_ε(ε_i)=σ^2_ασ_α>0
      η_i    α_i ϕ_i θ_i ε_i,tε_i.
   
  We further define 
  α̃_i=μ_α+δ_i'x_i, T_i^*+1. 
  We will investigate the post-shock aggregated estimators in M_2 
  in settings where δ_i is either fixed or random. 
  We let M_21 denote model M_2 with δ_i = δ for i= 1, …, n+1, 
  where δ is a  fixed unknown parameter.
  We let M_22 denote model M_2 with the following random effects 
  structure for δ_i:
  
      δ_iF_δ with E_F_δ(δ_i)=μ_δ, Var_F_δ(δ_i)=Σ_δ   with   δ_i  ε_i.
  
  We further define the parameter sets
  
      [                                                                         Θ                                                                         = {(η_i, ϕ_i, θ_i, α_i, x_i,t, y_i,t-1, δ_i) t= 1, …, T_i, i = 2, …, n +1},;                                                                       Θ_1                                                                         =          {(η_i, ϕ_i, θ_i, α_i, x_i,t, y_i,t-1, δ_i) t= 1, …, T_i, i = 1}, ]
  
  where Θ and Θ_1 can adapt to M_1 by dropping δ_i. We assume this for notational simplicity.
  
  
  
  
  
   §.§ Forecast
  
  
  In this section we show how post-shock aggregate estimators improve upon standard 
  forecasts that do not account for the shock effect.
  
  
  More formally, we will consider the following candidate forecasts: 
  
      Forecast 1: ŷ_1,T_1^*+1^1 = η̂_1 
          + ϕ̂_1 y_1,T_1^* + θ̂_1'_1,T_1^*+1 
          , 
         Forecast 2: ŷ_1,T_1^*+1^2 = η̂_1 
          + ϕ̂_1 y_1,T_1^* + θ̂_1'_1,T_1^*+1 
          + α̂,
  
  where η̂_1, ϕ̂_1, and θ̂_1 are all OLS estimators of η_1, ϕ_1, and θ_1, respectively, and α̂ is some form of estimator for the shock effect of time series of interest, i.e., α_1. 
  The first forecast ignores the presence of α_1 while the second forecast 
  incorporates an estimate of α_1 that is obtained from the other independent forecasts under study. 
  
  
  Note that the two forecasts do not differ in their predictions for 
  y_1,t, t = 1,… T_1^*. Instead, they only differ in predicting 
  y_1,T_1^*+1. Throughout the rest of this article we show that the donor pool of time series {y_i,t t = 1,…,T_i, i = 2,…,n+1} has 
  the potential to improve the forecasts for y_1, T_1^*+1 under  different 
  circumstances for the dynamic panel model M_1, M_21, and M_22. Improvement will be measured by assessing the reduction in risk that Forecast 2 offers over Forecast 1. We will return to the theoretical details of risk-reduction in Section <ref>.
  
  We specifically focus on predictions for y_1,T_1^*+1, the first post-shock response. It is important to note that in general α̂ does not converge to α_1 in any sense.  Despite this shortcoming, adjustment of the forecast for y_1,T_1^*+1 through the addition of α̂ has the potential to lower forecast risk under several conditions corresponding to different estimators of α_1. 
  
  
  
  
   §.§ Construction of shock effects estimators
  
  
  
  We now construct the aggregated estimators of the shock effects that appear in Forecast 2 (see Section <ref>). We use these to forecast response values y_1, T_1^*+1 assuming that T_1^* is known.  First, we introduce the procedures of parameter estimation for M_1, M_21, and M_22 (see Section <ref>). For i = 2, …, n+1, all parameters in this model will be estimated with ordinary least squares (OLS) using historical data of t = 1, …, T_i. For i = 1, we estimate all the parameters but α_1 using OLS procedures for t=1, …, T_1^*. In particular, let α̂_i, i = 2, …, n+1 be the OLS estimate of α_i.  Note that parameter estimation for M_1 is identically the same as that for M_21 or M_22.  We emphasize that α_is are random variables, but the OLS estimation is conditioned on the realizations from some distribution.
  
  
  Second, we introduce the candidate estimators for α_1. Define the 
  adjustment estimator for time series i=1 by
  
      α̂_adj = 1/n∑_i=2^n+1α̂_i,
  
  where the α̂_is in (<ref>) are OLS estimators of  
  the α_is for i = 2, …, n+1.  We can use α̂_adj as an estimator for 
  the unknown α_1 term for which no meaningful estimation information 
  otherwise exists. It is intuitive that α̂_ adj should perform 
  well under M_1 where we assume that α_i's share the same mean 
  for i= 1, …, n+1. However, it can also be shown that 
  α̂_ adj may be less favorable in M_21 
  and M_22, which will be discussed in detail in Section <ref>. 
  
  We also consider the inverse-variance weighted estimator 
  in practical settings where the T_i's and T_i^*'s vary greatly across i=2, …, n+1. 
  The inverse-variance weighted estimator is defined as 
  
      α̂_ IVW = ∑_i=2^n+1α̂_i / σ̂_iα^2/∑_i=2^n+1 1/σ̂_iα^2,
           where   σ̂_iα^2 = σ̂^2_i( 𝐔_i'U_i)_22^-1,
  
  where  α̂_i is the OLS estimator of α_i, 
  σ̂_i is the residual standard error from OLS estimation, 
  and U_i is the design matrix for OLS with respect to time series 
  for i = 2, …, n+1. Note that since σ is unknown, estimation 
  is required and the numerator and denominator terms are dependent in general. 
  
  
  
  However, α̂_ IVW can be a reasonable estimator in 
  practical settings. 
  
  
  
  We do not provide closed form expressions for α̂_ IVW 
  and α̂_ IVW but empirical performance of 
  α̂_ IVW is assessed via Monte Carlo simulation 
  (see Section <ref>).
  
  
  We now motivate a weighted-adjustment estimator for model M_21 
  and M_22. Our weighted-adjustment estimator is inspired by the 
  weighting techniques in synthetic control methodology (SCM) developed 
  in <cit.>. 
  
  
  However, our weighted-adjustment estimator is not a causal estimator and 
  our estimation premise is a reversal of that in SCM. 
  
  
  Our objective is in predicting a post-shock response y_1,T_1^*+1 that is not yet 
  observed using other time series whose post-shock responses are observed.
  
  
  
  We use similar notation as that in <cit.> to motivate our weighted-adjustment estimator. Consider a W∈^n weight vector 
  W=(w_2, …, w_n+1)', where w_i∈ [0,1] for all 
  i = 2, …, n+1. Construct
  
      X_1 = x_1, T_1^*+1',
          X = [     x_2, T_2^*+1';                 ⋮; x_n+1, T_n+1^*+1' ],
          and  X̂_1(W) 
          = W'X,
  
  where X_1,X̂_1(W) ∈^1 × p. Define W={W∈ [0,1]^n  1_n'W = 1 }. 
  Suppose there exists W^*∈W with 
  W^*=(w_2^*, …, w_n+1^*)' such that
  
      X_1=X̂_1(W^*),     i.e.,   x_1, T_1^*+1 = ∑_i=2^n+1 w_i^*x_i, T_i^*+1.
  
  Note that (<ref>) tries to find W^* such that x_1, T_1^*+1 is a convex combination of x_i, T_i^*+1 for i = 2, …, n+1 with weights W^*. Therefore, W^* should exist as long as X_1 falls in the convex hull of 
   
      {x_2, T_2^*+1', …, x_n+1, T_n+1^*+1' }.
  
  
  
  Our weighted-adjustment estimator will therefore perform well when the pool of time series posses similar covariates to the time series for which no post-shock responses are observed. We compute W^* as
  
      W^* = _W∈WX_1-X̂_1(W)_p.
  
  
  <cit.> commented that we can select W^* 
  so that (<ref>) holds approximately 
  
  and that weighted-adjustment estimation techniques of this form are not 
  appropriate when the fit is poor. 
  Note that W^* is not random since the covariates are assumed to be fixed. Since W is a closed and bounded subset of ^n,  W is compact. Because the objective function 
  is continuous in W, W^* will always exist. 
  
  Our weighted-adjustment estimator for the shock effect α_1 is
    
      α̂_ wadj = ∑_i=2^n+1 w_i^*α̂_i
             for   W^* = [   w^*_2       ⋯ w^*_n+1 ]'.
  
   
   We further define
  
      𝐕 = (𝐱_2, T_2^*+1, …,𝐱_n+1, T_n+1^*+1).
  
  
     If 𝐕 has full rank and it exists some 𝐖 satisfies (<ref>), the solution to  (<ref>) is unique.
   
  Proposition <ref> details some conditions when 𝐖^* is unique.  Note that 𝐕 is p × n. Therefore, if the covariates are of full rank and the true solution lies in the convex and compact 𝒲, a sufficient condition for 𝐖^* to be unique is p ≥ n. However, when p < n, 𝐖^* may not be unique. If it exists some 𝐖^* satisfies (<ref>) and p < n, there are infinitely many solutions to (<ref>).  The issue of non-uniqueness is further discussed in Section <ref>.
  
  
  
  In Section <ref> we specify that x_i,t, θ_i ∈^p. 
   
  However, it is not necessary that the all p covariates are important for every time series under study. The regression coefficients θ_i are nuisance parameters that are not of primary importance. 
  
  It will be understood that structural 0s in θ_i correspond to variables that are unimportant. 
  
  
  
    Our forecasting premise and estimation construction shares similarities with Bayesian viewpoints. From a Bayesian perspective, if we assign a prior π to α_1, α̂_ adj, α̂_ wadj, and α̂_ IVW can be interpreted as the Bayes rules with respect to π under different loss functions. If the sampling distribution of the data and π are known, it is possible to compute the Bayes risks of α̂_ adj, α̂_ wadj, and α̂_ IVW with respect to π,  thus enabling comparisons among them. Additionally, from Theorem 2.4 in Chapter 5 of <cit.>, α̂_ adj, α̂_ wadj, and α̂_ IVW are admissible if they are unique with probability one.
  
  
  
  
  § FORECAST RISK AND PROPERTIES OF SHOCK EFFECTS ESTIMATORS
  
  
  
  In this section, we discuss the properties that are related to forecast-risk reduction. In discussion of risk, it is useful to derive expressions for expectation and variance of the adjustment estimator α̂_ adj and weighted-adjustment estimator.  The expressions for the expectations are as follow,
  
   
      
    * Under M_1, α̂_ adj=α̂_ wadj = μ_α.
      
    * Under M_21, 
      α̂_ adj = μ_α + 1/n∑_i=2^n+1δ' x_i, T_i^*+1
      and
      α̂_ wadj = μ_α + δ'x_1, T_1^*+1 .
      
    * Under M_22,
      α̂_ adj = μ_α + 1/n∑_i=2^n+1μ_δ' x_i, T_i^*+1
      and
      α̂_ wadj = μ_α + μ_δ'x_1, T_1^*+1 .
    
  Formal justification for these results can be found in Appendix. Note that α̂_ adj, α̂_ wadj, and α̂_ IVW are not unbiased estimators for α_1. However, under M_1, α̂_ adj and α̂_ adj are unbiased estimators for α_1=μ_α. Nevertheless, α̂_ adj is a biased estimator for α_1 but α̂_ wadj is an unbiased estimator for α_1 under both M_21 and M_22. We collect these results in the following proposition. 
  
  
   
    
    
    * Under M_1, α̂_ adj is an unbiased estimator of α_1. Under M_21 and M_22, α̂_ adj is a biased estimator of α_1 in general.
    
    * Suppose that W^* satisfies (<ref>). Under M_1, M_21 and M_22, α̂_ wadj is an unbiased estimator of α_1.
  
  
  
  
  Unbiasedness properties for α_1 of α̂_ adj and α̂_ wadj allow for simple conditions for risk-reduction to hold, and more importantly motivates a bootstrap estimation for evaluation of these conditions. 
  These conditions and bootstrap will be discussed in Section <ref> and <ref>, respectively. Next, we present the variance expressions for α̂_ adj and α̂_ wadj as below.
  
  
    
    * Under M_1 and M_21,  
  
      α̂_ adj   =σ^2/n^2∑_i=2^n+1E{(U'_iU_i)^-1_22}+σ^2_α/n^2
      α̂_ wadj   = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22}+σ^2_α∑_i=2^n+1(w_i^*)^2
  
  
    * Under M_22, 
  
      α̂_ adj   =σ^2/n^2∑_i=2^n+1E{(U'_iU_i)^-1_22}+1/n^2(x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α)
      α̂_ wadj   = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22}
        + ∑_i=2^n+1 (w_i^*)^2 (x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α).
  
  
  Formal justification for these results can be found in Appendix. Note that the variances are not comparable in closed-form 
  because of the term E{(U'_iU_i)^-1_22}.  This term exists because of the inclusion of the random lagged response in our autoregressive model formulation.   
  
  
  
  Section <ref> details conditions needed for risk-reduction and comparisons of adjustment estimators. These conditions involve variances and expectations which may be difficult to compute in practice. To make use of those conditions in practice, estimation is required. Sections <ref> introduce a residual bootstrap procedure which estimates the involved parameters in those conditions and thus motivates prospective decision-making about whether α̂_i reduces the risk. Section <ref> describes our leave-one-out cross validation procedures, which prospectively estimate the correctness of such decision without observation of the post-shock response for the time series under study. Our simulations verify these procedures. 
  
  
  
  
   §.§ Risk-reduction conditions for shock effects estimators
  
  
  
  In this section we will discuss the conditions for risk reduction for individual shock effects estimators under M_1, M_21, and M_22. 
  For an adjustment estimator α̂, we will write the risk-reduction as Δ(α̂)=R_T_1^*+1,1-R_T_1^*+1,2 where R_T_1^*+1,2 is the risk of Forecast 2 calculated using the adjustment estimator α̂. 
  
  
  
    §.§.§ Conditions under M_1
  
   
   
  Recall that Proposition <ref> implies that the adjustment estimator α̂_ adj and weighted-adjustment estimator α̂_ wadj are unbiased for α_1 under M_1. With this result, we will have  the following propositions that specify the conditions that are necessary for risk reduction. 
  
  
  Under M_1,
  
    
    * Δ(α̂_ adj) > 0  when 
  (α̂_ adj) < μ_α^2.
    
    * if W^* satisfies (<ref>), then Δ(α̂_ wadj) > 0 when α̂_ wadj<μ_α^2. 
  
  
  
  
  Proposition <ref> says that under M_1 if the variance of the estimator is smaller than the squared mean of α_1, those estimators will enjoy the risk reduction properties.   In this setting, under M_1, Δ(α̂_ adj)= μ_α^2-α̂_ adj and Δ(α̂_ wadj)= μ_α^2-α̂_ wadj. 
  From Proposition <ref>, we obtain a risk-reduction condition
  
      α̂_ adj   =σ^2/n^2∑_i=2^n+1E{(U'_iU_i)^-1_22}+σ^2_α/n^2 < μ_α^2.
  
  
  Condition (<ref>) implies two facts: (1) adjustment (Forecast 2) is preferable to no adjustment (Forecast 1) asymptotically in n whenever μ_α≠ 0 (see Forecast in Section <ref>); (2) In finite donor pool settings, adjustment is preferable to no adjustment when μ_α is large relative to its variability and overall regression variability.   
  
  If  𝐖^* does not satisfy (<ref>), its unbiased properties for α_1 should hold approximately when the fit in (<ref>) is appropriate as commented in Section <ref>. From Proposition <ref> and the variance expression for α̂_ wadj, the risk-reduction condition for α̂_ wadj is
  
      α̂_ wadj
       = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22}+σ^2_α∑_i=2^n+1(w_i^*)^2 < μ_α^2.
  
  In this case, adjustment is preferable to no adjustment when μ_α is large relative to the weighted sum of variances for shock effects for other time series and overall regression variability.  However, the above criteria are generally difficult to evaluate in practice.  Sections <ref> and <ref>  provide detailed treatments on how to estimate the sign of Δ(α̂) in practice.
  
  
  
  
  
  
  
    §.§.§ Conditions under M_21 and M_22
  
  
  
  The shock effects α_is have different means under M_21 and M_22 unlike under M_1. 
  However, Proposition <ref> implies that α̂_ wadj is an unbiased estimator of α_1. 
  We now state conditions for risk-reduction.
  
  
   If W^* satisfies (<ref>), then Δ(α̂_ wadj) > 0 when α̂_ wadj<(α_1)^2  under M_21 and M_22. 
  
  Under Proposition <ref>, we can obtain a risk-reduction inequality that is similar to (<ref>),
  
      α̂_ wadj
       = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22} + ∑_i=2^n+1 (w_i^*)^2 (x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α) < (α_1)^2,
  
  where x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α may be replaced with σ^2_α in M_21. The conclusions and intuitions will be identically the same as what we have in Section <ref>. Proposition <ref> shows that α̂_ adj is a biased estimator of α_1 under M_21 and M_22 generally. Hence, Proposition <ref> no longer holds for α̂_ adj under M_21 and M_22.
  
  As an alternative, we can derive similar risk-reduction conditions that are appropriate for this setting. By Lemma <ref> (see Section <ref>) and risk decomposition, we will achieve risk-reduction as long as
  
      α_1^2= α_1+(α_1)^2
          >α̂_ adj-α_1^2
         =α̂_ adj +  (α̂_ adj-α_1)^2 
         =α̂_ adj +  α_1 + (α̂_ adj-α_1)^2.
  
  The above inequality simplifies to 
  
      (α_1)^2 >α̂_ adj+ (α̂_ adj-α_1)^2.
  
  
  
  
  
  
  
  
  
   
  
  
  As mentioned in Section <ref>, it is difficult to evaluate the expectation and variance of α̂_ IVW. We note that α̂_ IVW is generally biased for α_1. That is to say we can adapt the above proof to derive the risk-reduction conditions for α̂_ IVW: under M_1, M_21, and M_22, Δ(α̂_ IVW) > 0 when α̂_ IVW +(α̂_ IVW-α_1)^2<(α_1)^2.  In fact, more generally, using similar proof of Lemma <ref>, it can be shown that under M_2, the risk-reduction quantities are
  
      Δ(α̂_ adj) 
           = (α_1)^2 -α̂_ adj-(α̂_ adj-α_1)^2,
      Δ(α̂_ IVW) 
           = (α_1)^2 -α̂_ IVW-(α̂_ IVW-α_1)^2 ,
      Δ(α̂_ wadj) 
           = (α_1)^2 -α̂_ adj,
  
  where we estimate Δ(α̂) for estimator α̂  using bootstrap and leave-one-out cross validation procedures developed in Sections <ref> and <ref>. 
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
  
  
  
  
   §.§ Bootstrap for risk-reduction evaluation problems
  
  
  
  
  In this section, we present bootstrap procedures that approximate the distribution of our shock effect estimators, checks the underlying conditions of our risk reduction propositions, and estimate risk-reduction quantity using plug-in approach in practice. Our procedure involves the resampling of residuals in the separate OLS fits. This procedure has its origins in Section 6 of <cit.> and Chapter 12 of <cit.>. Our procedure involves the resampling of the residuals which are assumed to be the realizations of an iid process.
  
  Our first bootstrap procedure is as follows: Let B be the bootstrap sample size. At iteration b, first resample the indices I = {2, …, n+1} of the donor pool with replacement to form I^(b) with cardinality n, where we note that the elements of I^(b) may not be unique in terms of their  indices in the donor pool. Initialize y_i,0^(b)=y_i,0 for all i ∈ I^(b). Then, resample the residuals under models M_1 or M_2, compute the bootstrapped response  y_i,t^(b) for   t ∈{1, …, T_i} using the model estimated by original data,  and obtain shock effect estimators for each of the time series in the donor pool for all i ∈ I^(b). These shock effect estimators are then used to construct any of the adjustment estimators α̂^(b)_adj, α̂^(b)_wadj, and α̂^(b)_IVW, for b = 1,…,B. We can then estimate distributional quantities of our shock effect estimators under our considered models with the bootstrap samples α̂^(b)_adj, α̂^(b)_wadj, and α̂^(b)_IVW, for b = 1,…,B. We denote this procedure by B_u. We motivate a second bootstrap procedure B_f which treats the the donor pool as fixed, and not a realization from an infinite super-population. Therefore, there is no resampling of the donor pool in B_f, it is otherwise similar to B_u. An algorithmic formulation of B_u and B_f  are outlined in Section 2 in the Supplementary Materials.
  
  We will explicitly use these bootstrapped samples of shock effect estimators to check the risk-reduction conditions in Propositions <ref> and <ref>. Recall that α̂_ adj,  α̂_ wadj and α̂_ IVW are unbiased estimators of their expectations, and α̂_ wadj is an unbiased estimator of α_1 under M_1 and M_2 from Proposition <ref>. Our bootstrap procedure estimates the variance of our adjustment estimators. We can then estimate the risk-reduction propositions and inequalities. For example, we can estimate Δ(α̂_ adj) under model M_21 or M_22 with 
  
      Δ̂(α̂_ adj)     = (α̂_ wadj)^2 -S^2_α̂_ adj -(α̂_ adj-α̂_ wadj)^2,
  
  where S^2_α̂_ adj is the bootstrap sample variance estimator for (α̂_ adj). 
  
  
  
  
  
  
  
  
  We reiterate the philosophical distinction between B_u and B_f. B_u treats the donor pool as realizations from some infinite super-population of potential donors. In contrast, B_f treats the donor pool as being fixed  and known before the analysis is conducted, where the randomness arises from parameters and idiosyncratic error. A double bootstrap procedure with similar steps can estimate the distribution of Δ̂(α̂) for α̂∈A. The double bootstrap, instead of checking whether Δ(α̂)>0, can check whether a bootstrap percentile interval of resampled estimates of Δ(α̂) contain 0 at a desired error threshold. We investigated such a double bootstrap procedure and found that it produced inferences that were similar to those produced using the bootstrap techniques developed in the main text.
  
  We stress that our bootstrap approximations cannot alleviate the inherent bias of using our adjustment estimators as surrogates for α_1. We caution that the bootstrapping residuals in OLS estimation may not provide valid inference in moderate or high dimension where p < T_i but p / T_i is not close to zero for i∈{2, …, n+1} <cit.>; see alternatives for residual bootstrapping in linear models in <cit.>. 
  
  Recall that 𝐖^*  may not be unique if the conditions in Proposition <ref> are not satisfied. Non-uniqueness might be a concern theoretically. This is due to the fact that infinitely many different weights can lead to infinitely many non-unique α̂_ wadj's all targeting on the same α_1. For example, consider the case where the size of  donor pool to be 2, α̂_2=1, α̂_3=2, and there are two solutions to (<ref>), say, 𝐖^*_1=(1,0) and  𝐖^*_2=(0,1). In this scenario, the weighted adjustment estimator induced by 𝐖^*_1 has variance 1 whereas the one by 𝐖^*_2 has variance 2. Nevertheless, even if α̂_i has the same variance across i = 2, …, n+1, the same issue would occur if there were infinitely many W^* with different norms. It is possible to resolve this issue by selecting a unique weight W^* that optimizes a desirable objective function, prior to which one should find the bases spanning the subspace of W^* satisfying (<ref>). Simulations in the Section 4 of the Supplementary Materials provide some evidence that non-uniqueness of 𝐖^* is not problematic for inferences.
  
  
  
  
   §.§ Leave-one-out cross validation
  
  
  
  In this section, we adapt leave-one-out cross validation (LOOCV) to our estimation context in order to provide prospective evaluations of our adjustment techniques. Our proposed LOOCV procedure has its roots in Section 7.10 of <cit.>.  Recall in Section <ref> that we are given the data {(x_i,t, y_i,t)  i = 1, …, n+1, t = 1, …, T_i}, where {(x_1,t, y_1,t) t = 1, …, T_1} is the data of the time series of interest and the remaining observations form the donor pool. For iteration m ∈{1,…,n} of our LOOCV procedure, we set aside {(x_m + 1, t, y_m + 1, t)  t = 1, …, T_m+1} as the time series of interest, and construct a new donor pool {(x_i, t, y_i, t)  i ∈I_m, t = 1, …, T_i}, where I_m={2, …, n+1}∖{m+1}. Since the post-shock response y_m+1, T_m+1^*+1 is observed, we can evaluate the performance of our adjustment estimators and the original forecast made without adjustment (i.e., Forecast 1 in Section <ref>).
  
  LOOCV can be very computationally intensive when n is large, especially when combined with bootstrapping. To alleviate these concerns we can perform LOOCV with a random subset of k ≤ n iterations selected without replacement. In this setting, we let J be the randomly sampled indices. For m ∈J, we set aside {(x_m + 1, t, y_m + 1, t)  t = 1, …, T_m+1} as the time series of interest, and construct a new donor pool {(x_i, t, y_i, t)  i ∈I, t = 1, …, T_i}, where I={2, …, n+1}∖{m+1}. Based on the new donor pool, we estimate relevant parameters using bootstrap procedures outlined in Section <ref>. In other words, k times of bootstrapping are nested in a LOOCV procedure.  We find that k=5 or k=10 iterations of LOOCV performs well.
  
  We now outline how LOOCV can be used to prospectively assess the performance of adjustment estimators. Let A be the set of adjustment estimators. For each α̂∈A, let δ_α̂ = I(Δ̂(α̂)>0) be a decision rule where I(·) is the indicator function and a 1 corresponds to the decision to use  estimator α̂. If Δ(α̂)>0 (Δ(α̂)<0, respectively) but δ_α̂ incorrectly reported 1 (0, respectively) so that it makes the decision not to use α̂ (to use α̂, respectively), δ_α̂ is said to be incorrect. If Δ(α̂)<0 (Δ(α̂)>0, respectively) and  δ_α̂ correctly reported 0 (1, respectively) so that it makes the decision to use α̂, δ_α̂ is said to be correct. These situations are depicted in the following table: 
  
  
  
  
    
  
        
                  2cDecision 
  
                  δ_α̂ = 1     δ_α̂ = 0 
   
                  2*Truth      Δ(α̂)>0     Correct     Incorrect 
  
        3-4
            Δ(α̂)<0      Incorrect     Correct 
  
        
    
  
  
  
  
  
  
  
  
  We will use 
  C(δ_α̂) 
  	= I(δ_α̂ is correct)
  as a metric that evaluates the performance of forecasts made with the adjustment estimator α̂. If C(δ_α̂) > 0.5, we claim that δ_α̂ is better than random guessing. Note that C(δ_α̂) can generally be computed only when the post-shock response is observed. However, it is possible to estimate C(δ_α̂) using LOOCV. 
  The LOOCV estimates for C(δ_α̂) are
  
      C̅(δ_α̂)= 1/n∑_m = 1^n C^(-m)(δ_α̂),
  
  where C^(-m)(δ_α̂) is computed with respect to donor pool with index set I_m and the m+1 time series is treated as the time series of interest.  The LOOCV with k random draws estimates C(δ_α̂) as
  
      C̅^(k)(δ_α̂)= 1/k∑_m ∈JC ^(-m)(δ_α̂)	,
  
  where J is the set of the k randomly sampled indicies. 
  
  
    Note that we allow the time series within the donor pool to be dependent but donor pool should be independent of the time series of interest. However, if we assume the mutual  independence structure, C̅(δ_α̂) will be an almost unbiased estimator of C(δ_α̂) <cit.>. 
  
  
  
  
  
  
  § NUMERICAL EXAMPLES
  
  
  
  
  
   §.§ Modeling setup
  
  
  In this section we provide justification for our methods based on Monte Carlo simulation. We implemented our simulation based on M_22 with negligibly small Σ_δ approximating the design of ℳ_21.  We consider p=25 and μ_α=2, where p = 25 is set to satisfy conditions in Proposition <ref>. Parameter setup of our simulations is detailed as follows: the ϕ_i's are sampled independently from Uniform(0,1). We sampled T_i's independently from  Gamma(15, 10) that are further rounded to integers, where the minimum allowable value of T_i is fixed to be 90. We will randomly draw T_i^* from {p + 4, …, T_i-1}. The choices of T_i and T_i^* are set up to satisfy a necessary condition for the design matrix of OLS estimation to have full rank. Moreover, it is  designed to illustrate the performance of α̂_ IVW that may perform well in time series with varying lengths. Additionally, we generated the covariates from Gamma(1,2) to set up a setting when the α̂_ wadj may perform well. Last, we set δ_i10.5 and θ_i ∼01. We will consider parameter setup by varying σ in the model of y_i,t, n, the donor pool size, and σ_α in the model of α_i. We choose a Monte Carlo sample size of 30 replications and a bootstrap sample size of B = 200 for computation. Means and standard errors for estimated quantities will be recorded. Our LOOCV procedure will consider k=5 random draws. Recall in Section <ref> that B times of bootstrap are nested in a LOOCV with k random draws. It implies that B(k+1) times of bootstrap replications are required for each Monte Carlo simulation.
  
  
  
  
   §.§ Performance metrics
  
  
  Our adjustment estimators will be evaluated by multiple criteria. We interpret δ_α̂= I(Δ̂(α̂)>0) for α̂∈A as the guess, with 1 indicating that α̂ provides risk-reduction over the simple no-adjustment forecast, and 0 indicates the converse. We will consider the LOOCV estimators (<ref>) and (<ref>) to assess correct decision making. We will also consider the Euclidean distance between the post-shock forecasts ŷ_1,T_1^*+1, ŷ_1,T_1^*+1+ α̂_ adj, ŷ_1,T_1^*+1+α̂_ wadj, and ŷ_1,T_1^*+1+α̂_ IVW and the realized post-shock response y_1,T_1^*+1. The first two metrics can combine to assess our forecasting methodology prospectively while the latter requires the realization of the post-shock response y_1,T_1^*+1. 
  
  
  
  
  
   §.§ Monte Carlo results
  
  
  
  In this section, we discuss simulation results for the bootstrap procedures used in estimating parameters for risk-reduction propositions and inequalities. We mainly discuss simulations under M_2 (see Section <ref>) for B_u and B_f (see Section <ref>) with comparisons to those under M_1 whose results are listed in Section 3 in the Supplementary Materials. Two simulation setups are investigated. 
  
  In the first simulation setup, we consider the parameter combination of  n ∈{5, 10, 15, 25} and σ_α∈{5, 10, 25, 50, 100} where we fix σ=10. Note that α_1=52, where the last expectation is operated under the density of the covariates. In other words, data with σ_α∈{5, 10, 25, 50, 100}  should well represent the situations when the signal of the covariates is strong and when it is nearly lost. Results are displayed in Table <ref> in the Appendix, Section <ref>.
  
  
  In the second simulation setup, we consider the parameter combination of σ, σ_α∈{5, 10, 25, 50, 100} where we fix n=10. Likewise, σ, σ_α∈{5, 10, 25, 50, 100} will produce situations when the signal of the covariates is strong and when it is nearly lost in the model of  both y_i,t and α_i. Results are  in Table <ref> in the Appendix, Section <ref>.
  
  
  
  
  First, assuming that C̅^(k)(δ_α̂) well estimates C(δ_α̂) and fixing n, we  observe from Table <ref> that the decision making of δ_α̂ is nearly correct  for α̂∈A when σ_α is small  from Table <ref>. The reasons can be explained as follows. When σ_α is small, the signal of the covariates is strong so that α̂_ wadj will be expected to capture the signal according to construction of α̂_ wadj in Section <ref>. Moreover, when σ_α is small, M_22 approximates M_21 such that estimation of α_1 should be nearly unbiased according to Proposition <ref>. However, when the signal of the covariates is poor (σ_α is large), the decision rule δ_α̂  becomes unreliable for α̂∈A. It is to be expected since the bootstrap estimates become more biased. However, users can be warned by C̅^(k)(δ_α̂) to have an idea of the effectiveness of δ_α̂. Second, fixing σ_α, we can observe that the correctness of δ_α̂ increases when n increases. It is due to the robustness gain in estimation when n increases.  
  
  
  Additionally, we  observe that in most cases δ_α̂_ wadj reports α̂_ wadj reduces the risk even when C̅^(k)(δ_α̂_ wadj) starts to break down, though they  follow similar patterns.   Recall from Section <ref> that Δ̂(α̂) contains the squared bias for estimating α_1. But it is not present for Δ̂(α̂_ wadj) since we applied the fact  α̂_ wadj is unbiased for α_1 from Proposition <ref> in plugging it in with replacing α_1. Therefore, when the signal from covariates is poorer, δ_α̂_ wadj becomes less conservative. Besides, the averaged I(Δ̂(α̂)>0) times C̅^(k)(δ_α̂) can provide an approximation for the probability that α̂ actually reduces the risk assuming an symmetry of correctness between the cases when Δ̂(α̂)>0 and when Δ̂(α̂)<0. For example, when n = 5 and σ_α=50, the probability that α̂_ adj reduces the risk is approximately 0.83 × 0.59 = 0.490 from Table <ref>.  In other words,  the probability that α̂ reduces the risk has the same pattern as C̅^(k)(δ_α̂) has with n and σ_α for α̂∈A.
  
  From columns related to distance to y_1, T_1^*+1 in Table <ref>, as σ_α increases, the prediction appears to be poorer. When σ_α= 5, 10, 25, forecasts using α̂_ adj, α̂_ wadj, and α̂_ IVW are always better than the original forecast significantly. But it does not hold generally for the case when σ_α=50, 100. It is reasonable in that when the σ_α is large, it is difficult to find a reliable estimate of α_1.  Nevertheless, no statistical evidence has been found to support the claim that n matters in prediction. In other words, the size of the donor pool matters for producing reliable decision-making of δ_α̂ rather than reliable prediction. 
  
  
  
  
  
  
  From Table <ref>, we observe that as σ_α increases fixing σ, C̅(δ_α̂)   decreases, which is a pattern similar to the one shown in the first experiment. Furthermore, as σ increases fixing σ_α, C̅(δ_α̂) decreases as well. Note that the correctness  hinges on the estimation of the parameters. Since α̂_ wadj  is a linear combination of OLS estimates, as σ increases, α̂_ wadj increases as well. Therefore, α̂_ wadj become more volatile and its estimation of α_1 can be less reliable. Those reasons can explain  why an increase of σ_α contributes to a decrease of  C̅(δ_α̂).  We observe similar patterns for distance to y_1, T_1^*+1 as well. When σ increases with fixing σ_α, it is likely that the degree of variation of y_1,t exceeds the extent of adjustment improvement α̂ can contribute to for α̂∈A.
  
  
  
  
  With respect to averaged I(Δ̂(α̂)>0) (i.e., the guess), it starts to decrease as σ increases. This is reasonable if we believe the bootstrap estimate S^2_α̂ provides a good approximation for α̂ for α̂∈A. The reasons are outlined as follows: Recall in Section <ref>, the conditions of risk-reduction propositions involve (α_1)^2 > α̂ + (α̂-α_1)^2 for α̂∈A. Notice that α̂ is an increasing function of σ since α̂ is estimated by OLS. Therefore, it explains the reason why the increase of σ would result in a decrease of averaged I(Δ̂(α̂)>0) since the inequality is not likely to hold when  α̂ increases.
  
  
  Simulation for B_f with the same parameter setup as that of B_u are implemented. See Table <ref> and Table <ref> for results in the Appendix, Section <ref>. Comparing Table <ref>
   and Table <ref> yields that when n is moderately small (n = 10) and σ_α is small (σ_α=5), B_u is better than B_f with statistical evidence. For other situations, B_u  and B_f  are rather similar. It is likely that the extra randomness from sampling with replacement from donor pool compensates for the possible noises from a small donor pool. Concerning  Table <ref> and Table <ref>, it appears that when n = 10 and σ_α=5, B_u is better than  B_f  when σ increases. It might be the case that additional layer of bootstrap in the donor pool buffers the negative effects on C̅(δ_α̂)  introduced from increasing variation of y_i,t. However, when σ_α increases over 5 and n = 10,  B_f and B_u are quite similar under  situations of different σ and σ_α. In conclusion, B_u is better than B_f when the signal of the covariates is strong and n is moderately small; otherwise, they are similar. 
   
  
   Simulation results corresponding to model M_1 are listed in Section 3 in Supplementary Materials. Results under model M_1 are very similar to those of M_2, except for the difference among estimators. The results show that (1) the performance of α̂_ adj and  α̂_ IVW are nearly the same and (2) in many situations, α̂_ adj and  α̂_ IVW are  better than α̂_ wadj; in other situations, they are mostly the same. Recall that in M_1, the models for α_1 do not involve the covariates. Therefore, similarity weighting may not be informative when the model for α_i is identified wrongly. Under M_1, simple averaging, aimed for a reduction of variance, or inverse-variance weighting, targeting on reducing negative effects from varying time lengths, may work better.
  
  
   
  We have implicitly assumed that 𝐖^* is non-degenerate in the population in these simulations.  Recall that in Section <ref> we noted that if there exists some 𝐖^* which satisfies (<ref>) and p < n, then there will be infinitely many solutions to 𝐖^*. In this scenario, 𝐖^* will take values on the boundary of W, in which case bootstrapping may fail to estimate the distribution of α̂_wadj <cit.>.  When p < n and there exists some 𝐖^*∈W satisfies (<ref>), B_u fails since the non-uniqueness due to p < n will guarantee degeneracy of  𝐖^*. However, this issue will not occur under B_f since it takes W^* as being fixed and the parameter space is Θ that does not involve the constrained W. Moreover, simulations in Section 4 of the Supplementary Materials show that non-uniqueness does not seriously compromise the inference. 
  
  
  
  
  
  
  
  
  
  § FORECASTING CONOCO PHILLIPS STOCK IN THE PRESENCE OF SHOCKS
  
  
  
  We demonstrate our post-shock forecasting methodology on a time series of Conoco Phillips share prices after the occurrence of a structural shock. Conoco Phillips is a large oil and gas resources company <cit.>. The particular post-shock response that we predict happened after trading ended on Friday March 6th, 2020 and before trading began on Monday March 9th, 2020. It is reasonable that the timing of this shock is known, several events occurred over the trading weekend which had an impact on stock markets and the oil markets. For example, Russia and OPEC began a battle for global oil price control on Sunday, March 8th <cit.>, and several US states began declaring state of emergencies in response to the evolving coronavirus pandemic <cit.>. In this analysis we make the following design considerations:
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    (1) Selection of model. We will use an AR(1) model to forecast Conoco Phillips stock price. This model has been shown to beat no-change forecasts when predicting oil prices over time horizons of one and three months <cit.>.  For these reasons, we will consider 30 pre-shock trading days and we will forecast the immediate shock effect. All estimates will be adjusted for inflation. The model setup for AR(1) is exactly the same as what is stated in Section <ref> with addition of shock effects. All the parameters are estimated using OLS.
  
    (2) Selection of covariates. We consider different covariates for the model of α_i and y_i, t. The model of α_i incorporates daily S&P 500 index prices, West Texas Intermediate (WTI) crude oil prices, dollar index, 13 Week treasury bill rates, and Chicago Board Options Exchange volatility index (VIX). The model for y_i,t disregards VIX. This is because VIX is a metric for capturing market risk and sentiment, which is highly influential for the shock effect model.
  
    (3) Construction of donor pool. Our donor pool consists of Conoco Phillips shock effects observed in the past. We consider shock effects which occurred on March 14, 2008, several days in September, 2008, and November 27, 2014. The first sets of shock effects were observed during recessions that possessed similar characteristics to the current recession. In particular, all of these recessions were predicated by an inversion of the yield curve <cit.>. These 2008 shock effects correspond to the collapse of Bear Stearns, the placement of Fannie May and Freddie Mac in conservatorship on September 7th, the collapse of Lehman Brothers on September 15th, and the closing of Washington Mutual on September 25th <cit.>. The last shock effect corresponds to an OPEC induced supply side shock effect <cit.>. 
  
  
  
  
  
  
  
  
  We assume that the five shocks are independent of the shock that Conoco Phillips experienced on March 9, 2020. The covariates and response of time series in the donor pool are adjusted for inflation. Note that there are three shock effects nested in the time series 2008 September, we assume that these three shocks are independent, where the assumption checks using likelihood ratio test are provided in the Section 1 in the Supplementary Materials. The estimated shock-effects for α_i are -0.922, -7.063, -5.777, -6.395, -4.207 for i = 2, …, 6, respectively. Under M_2, we computed α̂_ adj, weighted adjustment α̂_ wadj, and α̂_ IVW. Note that non-uniqueness problems will not occur in this analysis since the conditions of Proposition <ref> are satisfied. To avoid the effect of  unit differences on weighting, we center and scaled the covariates in weights computation but not in the model of y_i,t. For α̂_ wadj, we observe that 
  𝐖^*= (0.000,0.000, 0.000, 0.273, 0.727)
  and 
  ‖X_1-X̂_1(W^*)‖_2 = 3.440.
  Note that the norm is computed using the k-dimensional Euclidean metric. The solution 𝐖^* suggests that the shock effect of interest is very similar to the September 25, 2008 shock effect and  the November 27, 2014 shock effect. 
  
  
  
  
  
  
   
   
   
   
  
  
  
  The resulting shock effect estimates are α̂_ adj=-4.872, α̂_ wadj= -4.805, and α̂_ IVW=-4.384. Using the bootstrap procedure B_f, we estimated parameters for risk-reduction propositions and risk-reduction quantities proposed in  Section <ref>. The estimated bootstrap variances for α̂_adj, α̂_wadj, and α̂_IVW  are 0.419, 0.559, and 0.667, respectively. 
  We verify the consistency of the result yielded by risk-reduction propositions with the reality as below.
  
  
  
  
  We can see from Figure <ref> that α̂_adj, α̂_wadj and α̂_IVW perform decently well. Our estimation framework misses the yet to be observed pot-shock response, however they perform much better than unadjusted forecasts that do not account for shock effects. The unadjusted forecast misses the post-shock response by 9.870 dollars whereas 
  the use of α̂_adj, α̂_wadj, and α̂_IVW misses by 5.324, 5.392, and 5.813 dollars, respectively. 
  
  The true shock effect is not fully recovered by α̂_adj, α̂_wadj, and α̂_IVW. This may be a result of a poorly constructed donor pool. The shock on March 9th, 2020 is in the midst of the COVID-19 pandemic and oil production volatility. It is difficult to find available stock market time series data that were generated under a similar setting. In any event, the shock on March 9th, 2020 was the largest price shock to Conoco Phillips shares by a wide margin, even after adjusting for inflation.
  
  
  From another perspective, it is possible that the stock of Conoco Phillips actually experienced multiple shocks on 2020 March 9th. For example, <cit.> studied the effect that different supply and demand shocks have on oil prices through a vector autoregressive model. Their model postulates an additive nature of shock effects, although the additivity parameters requires estimation in their context. Motivated by <cit.>, we also studied additive shock effect estimators where the shock effects corresponding to separate supply and demand shocks are added to estimate the unknown shock effect. The supply shock donor pool consists of the November 27, 2014 shock effect; and the demand shock donor pool consists of the remaining shock effects. The additive adjustment estimator computed by adding the α̂_adj, α̂_wadj, and α̂_IVW estimators for the demand and supply shock effects only, respectively, miss the post-shock Conoco Philips share price by 0.951, 0.405, and 1.460 dollars. These additive adjustment estimators do extremely well in this additive shock setting. 
  
  
  
  
  
  
  
  
  § DISCUSSION
  
  
  
  We developed a methodology for forecasting post-shock response values after the occurrence of a structural shock. Our methodology is as follows: construct a synthetic panel of time series which have undergone similar shocks, estimate the shock effects in those series, aggregate them, and then adjust the original forecast by adding the aggregated shock effect estimator to the original forecast. There have been several other recent similar methods developed for forecasting COVID-19 cases. For example, <cit.> constructed a Bayesian hierarchical model embracing data integration to improve predictive precision of COVID-19 infection trajectories for different countries. A similar setup may be appropriate for post-shock forecasting but may be too dependent upon model specification for the shock distribution. <cit.> employed a data-mining approach to combine COVID-19 data from different countries as input to predict global net daily infections and deaths of COVID-19 using a clustering approach. However, there is a tremendous amount of volatility in this form of COVID-19 data, and the fit of this prediction method may be improved with modeling structure or preprocessing of the donor pool. <cit.> proposed a model-free synthetic intervention method to predict unobserved potential outcomes after different interventions given a donor pool of observed outcomes with given interventions. They also provided useful guidelines for how to estimate the effects of potential interventions by giving recommendations for choosing the metric of interest, the intervention of interest, time horizons, and the donor pool. Although the methodology in <cit.> is quite general, there is no guarantee for theoretical properties in prediction without assuming any distributional structure. 
  
  
  
  We provided risk-reduction propositions and empirical tools that can prospectively assess the effectiveness of our adjustment strategies in additive shock effect settings. The model, under which we verify these claims, is a simple AR(1) model. Similar results can be obtained for more general models such as AR(p), vector autoregression, and generalized autoregressive conditional heteroskedasticity models. Generally, multiple shock effects can be nested within a time series; and time series in the donor pool can be dependent. As an example, we  considered a dependency structure for the September 2008 shock effects in our analysis of Conoco Phillips stock. But we note that consistency estimates from LOOCV may not work well if donor pool candidates are not mutually independent since the almost unbiased property hinges on the mutual independence among candidates in the donor pool.
  
  
  
  Although it is reflected in M_2, we stress that our proposed methods allow α_i to follow arbitrary distributions provided that its first and second moments exist. The covariates in the model for α_i under M_2 can be different from the covariates in the model of y_i,t.  Additionally, α_i can be heteroskedastic; i.e.,  they can have different variances. In this scenario, all the theoretical properties still hold though variance expressions in Section <ref> may not apply. We also note that our post-shock framework can be extended to settings where the shock effect can be decomposed into separable estimable parts. An example of this is the additive shock effect estimators that we studied in our Conoco Phillips analysis. Although our work is developed for time-series or AR(p) models, in fact, it can be generalized to any similar setting with a model of the response, whose parameters can be estimated unbiasedly, an additive shock-effect structure, and the structure that the time series in the donor pool are independent of the one of interest.
  
  
  
  
  Our bootstrap procedures can be extended to approximate the distribution of shock effect estimators from more general time series. If the data are subject to heteroskedasticity of unknown form, bootstrapping tuples of regressands and regressors proposed by <cit.> is robust in this situation with asymptotic validity in autoregressive models established by <cit.>. If serial correlation exists in the data, various block bootstrapping procedures <cit.> can be possible reasonable alternatives. Note that the pseudo time series generated by our proposed residual bootstrap are not stationary. If stationarity is of concern, one can be referred to the stationary bootstrap in <cit.> for stationary and weakly dependent time series. Nevertheless, it was shown that approximation accuracy might be a cost for the stationary bootstrap in autoregressive  models in finite sample <cit.>.  More work related to bootstrapping time series  can be referred to Chapters 3 and 4 in <cit.>, <cit.>, and Chapter 12 in <cit.>. It is up to  users in terms of selecting which procedure to choose but under different assumptions on the time series.
  
  
  
  Construction of donor pool is a critical step that directly matters in prediction. In our model, shocks from candidate time series in the donor pool, together with the shock in the time series of interest, should come from a family of distributions that allow for  varying variances and  varying means that are linear functions of common covariates. Remark <ref> shows that common covariates are not restrictive with structural zero in coefficients. However, it is difficult to verify such assumption or structure in practice. In principle, one should select identities that come from a common population, which is also the one  the identity under study come from, as candidates in donor pool. This selection criterion is often adopted in the literature of comparative case studies. For example, <cit.> used a set of cities in southern U.S. to estimate the effect of 1980 Mariel Boatlift on the Miami labor market. <cit.> utilized a donor pool of two Spanish regions to approximate the effect of terrorism on  the economic growth of Baseque Country. To estimate the effect of Proposition 99, a large tabacco control program in California  in 1988, on the annual per-capita cigarette sales of 2000, <cit.> constructed a donor pool consisting of US states, which did not implement any large-scale tobacco control program during their sampling period. Common features of those studies are that they all carefully constructed the donor pool and discarded any suspicious outliers. 
   
   
   
  
  Acknowledgements: We are grateful to Forrest W. Crawford, Karl Oskar Ekvall, Soheil Eshghi, Lutz Kilian, Ziyu Liu, Dootika Vats, and Dave Zhao for helpful comments.
  
  
  
  § APPENDIX
  
  
  
  
  
   §.§ Justification of Expectation of α̂_ adj and α̂_ wadj
  
  
  
  The building block for the following proof is the fact that least squares is conditionally unbiased conditioned on Θ. 
  
  Case I: under M_1: It follows that  under M_1 (see Section <ref>),
  
      α̂_ adj =1/n∑_i=2^n+1α̂_i|Θ = μ_α   and   α̂_ wadj   = ∑_i=2^n+1 w_i^* α̂_i|Θ= ∑_i=2^n+1 w_i^*μ_α= μ_α.
  
  where we used the fact that ∑_i=2^n+1 w^*_i=1. 
  
  Case II: under M_21 and M_22: Since ε̃_i, T_i=0, α̂̃̂_i=α̃_i=α_i, it follows that
    
      α̂_ wadj= E{E(∑_i=2^n+1 w_i^*α̂_i|Θ)}   = E( ∑_i=2^n+1 w_i^*α_i)
         = E{∑_i=2^n+1 w_i^*[μ_α+δ_i'x_i, T_i^*+1]}
         =μ_α+ μ_δ' ∑_i=2^n+1 w_i^*x_i, T_i^*+1W∈W
         =μ_α+  μ_δ'x_1, T_1^*+1. from (<ref>)
  
  Similarly,
    
      α̂_ adj   =μ_α+1/n∑_i=2^n+1μ_δ'  x_i, T_i^*+1.
  
  
  
  
  
   §.§ Justification of Variance of α̂_ adj and α̂_ wadj
  
  
  
  Notice that under the setting of OLS, the design matrix for M_2 is the same as the one for M_1. Therefore, it follows that
    
      α̂_ wadj   = α̂_ wadj|Θ + α̂_ wadj|Θ
         =E{Var(∑_i=2^n+1 w_i^*α̂_i|Θ)} +Var(∑_i=2^n+1 w_i^*α_i)
  
  Under M_21 where δ_i=δ are fixed unknown parameters,  we will have
    
      α̂_ wadj   = E{∑_i=2^n+1(w_i^*)^2(σ^2(U'_iU_i)^-1_22)} +σ^2_α∑_i=2^n+1(w_i^*)^2  
         = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22}+σ^2_α∑_i=2^n+1(w_i^*)^2.
  
  Similarly, under M_22 where we assume δ_i ε_i,t, we have
   
      α̂_ wadj   = σ^2∑_i=2^n+1(w_i^*)^2E{(U'_iU_i)^-1_22}
        + ∑_i=2^n+1 (w_i^*)^2 (x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α)
  
  For the adjustment estimator, we simply replace W^* with 1/n1_n. Thus, under M_21 we have 
   
      α̂_ adj   =σ^2/n^2∑_i=2^n+1E{(U'_iU_i)^-1_22}+σ^2_α/n^2
  
  Under M_22, we shall have
   
      α̂_ adj   =σ^2/n^2∑_i=2^n+1E{(U'_iU_i)^-1_22}+1/n^2(x_i, T_i^*+1'Σ_δx_i, T_i^*+1 + σ^2_α).
  
  Notice that M_1 differs from M_21 only by its mean parameterization of α (see Section <ref>). In other words, the variances of α̂_ adj and α̂_ wadj under M_1 are the same for those under M_21.
  
  
  
   §.§ Proofs for lemmas and propositions
  
  
  [<ref>]
    The proof of <cit.> in Appendix A.2 and A.3 adapts easily to Proposition <ref>.
  
  
  
  [<ref>] The proof for unbiasedness follows immediately from discussions related to expectation in Section <ref>. For the biasedness of  α̂_ adj under M_21 and M_22, we write the bias term for α̂_ adj as below.
  
      Bias(α̂_ adj) = 
        1/n∑_i=2^n+1δ'(x_i, T_i^*+1-nx_1, T_1^*+1)      for M_21
      1/n∑_i=2^n+1μ_δ'(x_i, T_i^*+1-nx_1, T_1^*+1)      for M_22.
  
  But it may be unbiased in some special circumstances when the above bias turns out to be 0. 
  
  
     The forecast risk reduction is R_T_1^*+1,1-R_T_1^*+1,2=α_1^2-α̂-α_1^2 for all estimators of α_1 that are independent of Θ_1 (see Section <ref>).
  
  
    
  [<ref>]
    Define 
    
      C(Θ_1) =η̂_1 +ϕ̂_1 y_1, T_1^*+θ̂_1'x_1, T_1^*+1 -(η_1 +ϕ_1y_1,T_1^*+θ_1'x_1,T_1^*+1),
  
    where Θ_1 is as defined in (<ref>). Notice that
    
      R_T_1^*+1,1= E{(C(Θ_1)-α_1)^2}     and       R_T_1^*+1,2= E{(C(Θ_1)+α̂-α_1)^2}.
  
    It follows that
    
      R_T_1^*+1,1-R_T_1^*+1,2=α_1^2-2C(Θ_1)α̂-α̂-α_1^2.
  
    Assuming S=(1_n, y_1,t-1, x_1) has full rank, under OLS setting, η̂_1, ϕ̂_1, and θ̂_1 are unbiased estimators of η_1, ϕ_1, and θ_1, respectively under conditioning of Θ_1. Since we assume α̂ is independent of Θ_1, through the method of iterated expectation,
    
      C(Θ_1)α̂=E{α̂·C(Θ_1)|Θ_1}=0.
  
    It follows that
    
      R_T_1^*+1,1-R_T_1^*+1,2=α_1^2-α̂-α_1^2,
  
    which finishes the proof.
  
  
  
  
  
  
  [<ref>] The proofs are arranged into two separate parts as below.
  
   Proof for statement (i): Under M_1, α̂_ adj is an unbiased estimator of α_1 because
    
      E( 1/n∑_i=2^n+1α̂_i)
         = 1/n∑_i=2^n+1α̂_i   = 1/n∑_i=2^n+1α̂_i|Θ
         =  1/n∑_i=2^n+1α_i
         = μ_α=α_1,
  
    where we used the fact that OLS estimator is unbiased when the design matrix U_i is of full rank for all i = 2, …, n+1. Because α_1ε_i,t, α̂_ adjα_1=α̂_ adjα_1=(α̂_ adj)^2. By Lemma <ref>, 
      
      R_T_1^*+1,1-R_T_1^*+1,2   =α_1^2-α̂_ adj-α_1^2
          =α_1^2-α_1^2- α̂_ adj^2+2α̂_ adjα_1
         = μ_α^2 - α̂_ adj
  
    Therefore, as long as we have α̂_ adj<μ_α^2, we will achieve the risk reduction. 
  
   Proof for statement (ii): By Proposition <ref>, the property that α̂_ wadj is an unbiased estimator of μ_α holds for M_1. The remainder of the proof follows a similar argument to the proof of statement (i).
  
  
  
  
  [<ref>]
    By Proposition <ref>, the property that α̂_ wadj is an unbiased estimator of α_1 holds for M_21 and M_22. The remainder of the proof follows a similar argument to the proof of Proposition <ref>.
  
  
  
  
  
  
  
   §.§ Tables for Section <ref>
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  plainnat
  
  
  
  
  
  