Shrinkage Estimation and Forecasting in Dynamic Regression Models Under Structural Instability

1 Introduction

A sizeable strand of the literature in economics and finance focuses on autoregressive (AR) models. These models are heavily used in forecasting economic and financial variables and are frequently considered as benchmarks in forecast competitions, as they are difficult to beat. Nevertheless, since many time series data in economics and finance are characterized by parameter instability, which is now widely recognized as an important source of forecast failure as recorded by Stock and Watson (1996), Hansen (2001), Giacomini and Rossi (2009), Rossi and Sekhposyan (2010), Inoue and Rossi (2011), Clements and Hendry (2006, 2011, and Rossi (2013) inter alia, there is an increasing evidence that parameters of AR models in many economic and financial time series are unstable and subject to structural breaks. For example, Mankiw and Miron (1986) and Mankiw, Miron, and Weil (1987) considered AR(1) models and find parameter instability in the short-term interest rate. Phillips, Wu, and Yu (2011) find parameter instability for 1990s NASDAQ stock prices. See also Garcia and Perron (1996) and Stock and Watson (1996) who document instability related to the autoregressive terms in a wide variety of economic time series. This suggests a need to study the forecasting performance of AR models when they are subject to structural breaks.

This paper considers an ARX model that contains structural breaks, and proposes a Stein-like shrinkage estimator that exploits observations in neighboring regimes.^[1] Our proposed estimator is a weighted average of a restricted estimator (which uses full-sample of observations under the restriction of no breaks in the parameters) with an unrestricted estimator (which uses observations within each regime). The restricted estimator is inconsistent when there is a break while it is the most efficient. On the other hand, the unrestricted estimator is consistent but the consistency comes at the cost of losing efficiency. Hence, the proposed shrinkage estimator trades off an increased bias introduced by the restricted estimator against a reduction in error variance resulting from using a full-sample of observations. The averaging weight depends on the weighted distance of the restricted and unrestricted estimators, which is similar to the James-Stein weight, cf. Stein (1956) and James and Stein (1961).^[2] Therefore, it assigns appropriate weights to each of the two estimators by measuring the magnitude of the structural break. We derive the analytical large-sample approximation of the bias, mean squared error (MSE) and risk for our proposed shrinkage estimator, the unrestricted estimator, and the restricted estimator.^[3] We derive the condition under which the risk of the shrinkage estimator is lower than the risk of the unrestricted estimator under any break size and break points. We also show how the method can be used in out-of-sample forecasting. Furthermore, we extend the results to the model with unit root and non-stationary regressors.

To the best of our knowledge, this is the first paper that analytically derives the asymptotic approximation of the bias, and risk up to order T ⁻¹, and MSE up to order T ⁻² in ARX models under structural breaks, where T is the total number of observations. Furthermore, we provide new results to the model with unit root and non-stationary regressors under structural breaks, which have not been considered before in the literature. Because of these considerations, we would like to point out that our results differ from the recent work by Lee, Parsaeian, and Ullah (2022a, 2022b, 2022c) who consider different types of combined estimators to improve forecasts under structural breaks in a static stationary time series model without lagged dependent variable and using different weights. For example, in Lee, Parsaeian, and Ullah (2022a) the combination weight is different and set to be a constant between zero and one, Lee, Parsaeian, and Ullah (2022b) consider Stein-like combined estimator, however, the optimal weights in this paper are different due to the presence of lagged dependent variable and theoretical frameworks. Lee, Parsaeian, and Ullah (2022c) consider another type of combined estimators which assigns a full weight of one to the post-break sample observations and a weight between zero and one to the pre-break sample observations. In addition, there are three additional main differences that have not been discussed in the previous works by Lee, Parsaeian, and Ullah (2022a, 2022b, 2022c). First, this paper considers a dynamic ARX model and allows for unit roots (or integrated of higher order) and non-stationary regressors. Second, in this paper the dominance property of the proposed Stein-like shrinkage estimator holds for any fixed deviations from the restrictions. This complements the “local asymptotic” argument discussed in Lee, Parsaeian, and Ullah (2022b). Third, we derive the analytical large sample approximation, using Nagar (1959) method, of the bias and risk up to order T ⁻¹, and MSE up to order T ⁻² for the proposed shrinkage estimator, the unrestricted estimator, and the restricted estimator under structural breaks. This allows us to theoretically and numerically compare the performance of these estimators. Because of the presence of the lagged dependent variable in the model considered in this paper, the theoretical results presented here are noticeably different than those in Lee, Parsaeian, and Ullah (2022a, 2022b, 2022c). Even in the special case of no lagged variable, our results are different with them because the combined estimator in Lee, Parsaeian, and Ullah (2022a) has different weights, and in Lee, Parsaeian, and Ullah (2022b, 2022c) they either derive only local asymptotic theory results instead of large sample approximation results and/or their weights in combined estimator are different.

We present an extensive Monte Carlo simulation study to evaluate the finite sample forecasting performance of the proposed shrinkage method. The results support our theoretical findings, and show the outperformance of our shrinkage estimator relative to the unrestricted estimator in finite sample. In particular, our numerical results show the benefits of exploiting the neighboring observations in estimation and forecasting relative to using only the observations within each regime. Furthermore, we undertake an empirical analysis for forecasting the output growth using 131 macroeconomic and financial time series to compare the forecasting performance of the shrinkage estimator with the unrestricted estimator, the restricted estimator, and a range of alternative methods existing in the literature. The empirical results suggest the out-performance of our method relative to the alternative methods, in the sense of mean squared forecast errors (MSFE).

Analysis of AR models subject to structural breaks have been also considered by Clements and Hendry (1998, 1999 and Pesaran and Timmermann (2005), but they focus on the analysis of forecast errors decomposition, while the focus of this paper is developing an estimator to deal with the parameter instability. Specifically, Clements and Hendry (1998, 1999 analyze the forecast errors from AR models subject to structural change by assuming that the parameters of the AR models remain unchanged during the estimation period. Pesaran and Timmermann (2005) consider the small sample properties of forecasts from AR models estimated from windows of different sizes. See also Banerjee and Urga (2006) for a comprehensive review of developments in the fields of modelling structural breaks.

The analysis of approximating the moments of estimators in autoregressive models dates back to Bartlett (1946), who finds a first-order variance approximation in an autoregressive gaussian process (see also Hurwicz 1950). White (1961) and Shenton and Johnson (1965) give approximations of the first two moments in the AR(1) model. Kendall (1954) and Marriott and Pope (1954) consider AR(1) models with intercepts, and find the approximate bias of the least-squares estimator of the lagged dependent variable coefficient. Recently, a number of papers have studied the small sample bias of the ordinary least-squares estimator in single dynamic regression models, see for example Grubb and Symons (1987) and Kiviet and Phillips (1993). Kiviet and Phillips (1994) extend the analysis of Kiviet and Phillips (1993) to the higher-order dynamic regression models. More recently, Kiviet and Phillips (2012) found the higher-order approximate bias of the least-squares estimator of the slope coefficients in a stable ARX(1) model. Kiviet and Phillips (2005) extend these results to non-stable ARX(1) models, and examine the moments of the least-squares estimator in the single normal ARX(1) model with an arbitrary number of exogenous regressors when the true coefficient of the lagged-dependent variable is unity.

The paper is organized as follows. Section 2 describes the model. For simplicity, we discuss the problem under a single break, which simplifies the essential idea without complicating notation. However, the generalization of the method for the multiple breaks is straightforward. In Section 3, we introduce the estimators. We give the bias, MSE, and the risk of the Stein-like shrinkage estimator using large-sample approximations in Section 4, while those of the unrestricted estimator are provided in the Supplementary Online Appendix B. We extend the results of Section 4 to a model where the slope coefficient of the lagged dependent variable is unity in Section 5. Monte Carlo results are given in Section 6. Results from our empirical example are given in Section 7. Conclusions are given is Section 8. Proofs and detailed calculations are provided in Appendix A.

Notation: Throughout the paper we adopt the following notation. I _p and 0 _p×q denote the p × p identity matrix and p × q matrix of zeros, respectively. tr(⋅) denotes trace, ⊗ denotes Kronecker products, and 1 ( ⋅ ) denotes the indicator function. For an m × n real matrix A = (a _ij ), we write the transpose A′. We write A = O(b) when the non-zero elements of A are of order b, i.e. a _ij = O(b). For a stochastic matrix A, we write A = O _p (b) when the non-zero elements of A are of order b, i.e. a _ij = O _p (b).

2 The Model

Consider the following first-order linear dynamic regression model (ARX(1) model)^[4] defined over the period t = 1, 2, …, T, which is subject to a single structural break at time T ₁

where y _t is the dependent variable, x _t = (x _t,1, …, x _t,k )′ is a k × 1 vector of exogenous regressors, and u _t is the unobserved error term.^[5] The (k + 1) × 1 vectors of the pre-break and post-break slope coefficients are denoted by α 1 = λ 1 , β 1 ′ ′ , and α 2 = λ 2 , β 2 ′ ′ , respectively, which are the parameters of interest. In a matrix form, the model can be expressed as

where y = (y ⁽¹⁾′, y ⁽²⁾′)′ is a T × 1 vector of the dependent variables, y ( 1 ) = ( y 1 , … , y T 1 ) ′ , y ( 2 ) = ( y T 1 + 1 , … , y T ) ′ . The T × 2(k + 1) matrix of observations on the regressors is denoted by Z = diag(Z ₁, Z ₂), where Z i = ( y − 1 ( i ) , X i ) for i = 1, 2, y − 1 ( 1 ) = ( y 0 , … , y T 1 − 1 ) ′ , y − 1 ( 2 ) = ( y T 1 , … , y T − 1 ) ′ , X 1 = ( x 1 , … , x T 1 ) ′ , and X 2 = ( x T 1 + 1 , … , x T ) ′ . u = (u ₁, …, u _T )′ is a T × 1 vector of disturbances, and α = α 1 ′ , α 2 ′ ′ is a 2(k + 1) × 1 vector of the unknown slope coefficients.

We note that u ₀ is the start-up error term, i.e. the error term at time 0. d = 0 represents the fixed start-up, and if d ≠ 0 the start-up is random. Assumption 1(ii) excludes non-stationary regressors including deterministic or stochastic trends, but the presence of such variables will not change the approximation bias, MSE, MSFE, and risk formulas in Section 4. Their inclusion only reduces the order of magnitude of the moments and the order of their remainder terms (for a similar discussion see Kiviet and Phillips 2012). We demonstrate in Section 5 that our results are still applicable under a unit root (relaxing Assumption 1(i)) and non-stationary regressors (relaxing Assumption 1(ii)). Assumption 1(v) assumes that the error terms are normally distributed and excludes serial correlation in the errors. This assumption can be relaxed at the expense of extra terms to be added in the bias and MSE of the estimators. However, since the model in (2.1) has a dynamic structure, it captures much of the serial correlation.

We analyze the model conditional on the observed matrix X in Section 4. Therefore, in order to distinguish the fixed and zero mean stochastic elements of the regressor matrix Z _i , i = 1, 2, we decompose Z i = Z ̄ i + Z ̃ i , where Z ̄ i is defined as the expectation of Z _i conditional on X and y ̄ 0 . For notational simplicity, we denote E ( ⋅ ) ≡ E ( ⋅ | X , y ̄ 0 ) . Hence,

and

where ς ₁ = (1, 0, …, 0)′ is a vector of (k + 1) × 1, and T ₂ = T − T ₁ is the number of post-break sample observations. Define the T _i  × T matrix Λ _i = (F _i , Δ _i ) for i = 1, 2. Then, we define y ̄ − 1 ( i ) = Λ i y ̄ − 1 * , where y ̄ − 1 * = y ̄ 0 , x 1 ′ β 1 , … , x T 1 ′ β 1 , x T 1 + 1 ′ β 2 , … , x T − 1 ′ β 2 ′ , F 1 = 1 , λ 1 , … , λ 1 T 1 − 1 ′ is T ₁ × 1, F 2 = λ 1 T 1 , λ 1 T 1 λ 2 , … , λ 1 T 1 λ 2 T 2 − 1 ′ is T ₂ × 1, and

are T ₁ × (T − 1) and T ₂ × (T − 1), respectively. Moreover, define the (T + 1) × 1 random vector ν = u 0 , u ′ ′ = ( u 0 , u 1 , … , u T ) ′ ∼  N ( 0 , Ω ν ) , where Ω ν = diag σ 1 2 I T 1 + 1 , σ 2 2 I T 2 , and the T × (T + 1) matrix G = G 1 ′ , G 2 ′ ′ , with the T _i  × (T + 1) matrix G _i = (dF _i , C _i ), and C i = ( Δ i , 0 T i × 1 ) , for i = 1, 2. Then, it can be easily verified that Z ̃ i = G i ν ς 1 ′ for i = 1, 2. Therefore, we have Z = Z ̄ + Z ̃ , where Z ̄ = diag ( Z ̄ 1 , Z ̄ 2 ) , and Z ̃ = diag ( Z ̃ 1 , Z ̃ 2 ) = ∑ i = 1 2 L i G ν e 1 , i ′ , where the 2(k + 1) × 1 vectors e 1,1 = ( ς 1 ′ , 0 1 × ( k + 1 ) ) ′ and e 1,2 = ( 0 1 × ( k + 1 ) , ς 1 ′ ) ′ , and L 1 = I T 1 0 T 1 × T 2 0 T 2 × T 1 0 T 2 × T 2 , L 2 = 0 T 1 × T 1 0 T 1 × T 2 0 T 2 × T 1 I T 2 , are T × T selection matrices.

3 Estimation

Our goal is to estimate the vector of slope parameters, α, in Equation (2.2). We consider three estimators of the slope parameters: (i) an unrestricted least squares estimator that estimates the slope parameters of each regime only using the observations within the same regime, (ii) a restricted estimator that shrinks the unrestricted estimator towards a restricted parameters space which ignores the break in the slope coefficients, and (iii) a Stein-like shrinkage estimator which is a weighted average of the restricted estimator and the unrestricted estimator, and the averaging weight is proportional to a weighted quadratic loss function.

4 Large-Sample Approximation of the Shrinkage Estimator

We employ the large-sample approximation method developed by Nagar (1959), to analyze the bias, MSE, and risks of the shrinkage estimator (conditional on X) under Assumption 1. To find moment approximations using the Nagar approach, we begin by expressing the estimation error in term of stochastic components which are of decreasing order of magnitude in terms of the sample size.^[9]

Theorem 1 gives the bias, MSE, and Risk of the Stein-like shrinkage estimator. The bias and MSE of the unrestricted estimator is given in Lemma B.3 in the Supplementary Online Appendix. We note that the risk of the Stein-like shrinkage estimator is generalized to allow for any positive definite weight matrix, W, of order O(T). Two arbitrary choices of W are TI _2(k+1), and W _f defined in Section 4.1, where the former provides an unweighted mean squared error, and the latter gives the mean squared forecast error studied in the next section.

In the following Corollary we show the dominance conditions of the Stein-like shrinkage estimator relative to the unrestricted estimator. Let ϕ W = α ′ P ̄ ′ W P ̄ α , μ = tr ( W P ̄ Q Σ ) + α ′ P ̄ ′ W Θ + tr ( W Ψ ) , and η = 2 α ′ P ̄ ′ Q − 1 P ̄ Q Σ W P ̄ α + 2 α ′ P ̄ ′ Q − 1 Ψ W P ̄ α + α ′ P ̄ ′ Φ W P ̄ α .

Corollary 1.1 shows that the proposed shrinkage estimator dominates the unrestricted estimator in terms of having a smaller risk when the shrinkage parameter satisfies the condition 0 < τ ≤ 2 ϕ ϕ W μ − 1 ϕ η . In addition, as the choice of the shrinkage parameter is user-specified, its optimal value and our ideal choice that minimizes the risk of the Stein-like shrinkage estimator up to O(T ⁻¹), is τ _opt. This implies that the Stein-like shrinkage estimator always performs better or performs equal to the unrestricted estimator which is one of the common methods of estimating the coefficients under structural breaks.^[12]

Since τ _opt depends on σ 1 2 , σ 2 2 , Q , Z ̄ , C , G , ω , α , and y ̄ 0 , which are unobserved, we consider the estimated optimal shrinkage parameter, denoted by τ ̂ opt , defined as

In Equation (4.7), μ ̂ , η ̂ , ϕ ̂ and ϕ ̂ W correspond to μ, η, ϕ and ϕ _W respectively, after replacing σ i 2 , α , Q , Z ̄ , C , G with their estimates, denoted by σ ̂ i 2 , α ̂ , Q ̂ , Z ̂ , C ̂ , G ̂ , where α ̂ and σ ̂ i 2 are the unrestricted estimators, Q ̂ = ( Z ′ Z ) − 1 , Z ̂ = diag ( Z ̂ 1 , Z ̂ 2 ) , with Z ̂ i = ( F ̂ i y 0 + C ̂ i X i β ̂ i , X i ) , C ̂ and F ̂ correspond to C = C 1 ′ , C 2 ′ ′ and F = F 1 ′ , F 2 ′ ′ after replacing λ ₁, and λ ₂ by λ ̂ 1 , and λ ̂ 2 . Similarly, Θ ̂ , Ψ ̂ , and P ̂ correspond to Θ, Ψ, and P ̄ after replacing the unobserved parameters with their estimates. Since we define the Stein-like shrinkage estimator with positive shrinkage parameter in (3.5), we set τ ̂ opt to zero when the condition μ ̂ > 1 ϕ ̂ η ̂ does not hold. In this case, the Stein-like shrinkage estimator assigns a weight one to the unrestricted estimator and a zero weight to the restricted estimator.

In the following corollary, we show that the estimated optimal shrinkage parameter, τ ̂ opt , is an unbiased estimator of the infeasible optimal shrinkage parameter τ _opt up to order T ⁻¹. Hence, when the sample size is large enough the risk of the Stein-like shrinkage estimator using the estimated optimal shrinkage parameter is smaller than the risk of the unrestricted estimator.

5 Unit Root

In this section, we extend the results of Theorem 1 by relaxing Assumption 1(i) and (ii) on the stability, and the stationarity of the exogenous regressors. We assume that for j = 1, …, k + 1, the series of real positive constants g _j is given, such that Z ̇ ′ Z ̇ = O p ( T ) , where Z ̇ = Z N , and N = I 2 ⊗ diag T − g 1 , … , T − g k + 1 .