10 带时间序列误差的回归模型

在统计学的数据分析中, 线性回归分析是最常用的分析工具之一。 线性回归以一元线性回归为例,模型如下

\[\begin{align} Y_t = \beta_0 + \beta_1 X_t + e_t, \ t=1,2,\dots,T \tag{10.1} \end{align}\]

其中自变量\(\{ X_t \}\)为常数列, \(\beta_0, \beta_1\)为未知的系数, \(\{ e_t \}\)为零均值独立同分布随机误差序列, 方差为\(\sigma_e^2\), 因变量\(\{ Y_t \}\)为随机变量列。 参数\(\beta_0, \beta_1, \sigma_e^2\)可以用最小二乘法估计, 估计量无偏、相合, 系数的估计渐近正态分布。 为了对估计结果做假设检验或者用模型结果做预测, 一般还假设\(\{ e_t \}\)为零均值独立同正态分布序列。

在金融研究中, \(\{ X_t \}\)\(\{ Y_t \}\)一般都是时间序列, 而且\(\{ e_t \}\)也是时间序列,有序列相关性。 这时, 最小二乘估计不是最有效的估计甚至于可能不相合, 基于\(\{ e_t \}\)独立同正态分布所做的标准误差估计、假设检验和预测都不再成立。

(Cochrane and Orcutt 1949)的文章指出, 当\(\{ e_t \}\)彼此正相关时, 回归系数的标准误差估计偏低, 使得相应的t和F检验统计量的绝对值偏大。 (Durbin and Watson 1950)(Durbin and Watson 1951)给出了一阶自相关的检验。 Durbin-Watson序列自相关检验使用检验统计量 \[ \text{DW} = \frac{\sum_{t=2}^T(\hat e_t - \hat e_{t-1})^2}{\sum_{t=1}^T \hat e_t^2} \approx 2(1 - \hat\rho_1) \] 其中\(\hat e_t\)是回归残差, 当无序列相关时DW统计量接近于2。 此统计量只用到一阶自相关。

\(\{ e_t \}\)是平稳可逆ARMA序列时, 可以将线性回归模型与平稳可逆ARMA序列同时估计, 可以得到需要标准误差估计、假设检验和预测。 R的arima()函数提供了一个xreg=用来引入回归自变量。

如果不关心\(\{ e_t \}\)的具体模型, 而只关心对回归系数的SE的正确估计以及假设检验的正确性, 可以仅假设\(\{ e_t \}\)的协方差结构而不考虑\(\{ e_t \}\)的建模。 这样的方法在混合线性模型中使用。

下面举例说明。 考虑美国国债一年期与三年期利率, 数据为周数据,从1962-01-05到2009-10-04,共2467个观测。 记一年期利率序列为\(\{ x_{1t} \}\), 三年期利率序列为\(\{ x_{3t} \}\)

先读入数据:

da <- read_table2(
  "w-gs1yr.txt",
  col_types=cols(.default=col_double()))
x1 <- da[["rate"]]
y1 <- xts(x1, make_date(da[["year"]], da[["mon"]], da[["day"]]))
da <- read_table2(
  "w-gs3yr.txt",
  col_types=cols(.default=col_double()))
x3 <- da[["rate"]]
y3 <- xts(x3, make_date(da[["year"]], da[["mon"]], da[["day"]]))
rm(da)
xts.gs <- cbind(gs1yr=y1, gs3yr=y3)
rm(y1, y3)

两个利率序列的时间序列图:

plot(xts.gs, type="l", 
     main="US Federal Bonds 1 Year and 3 Years",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto",
     col=c("blue", "red"))
美国一年期与三年期国债利率周数据

图10.1: 美国一年期与三年期国债利率周数据

蓝色为一年期,红色为三年期。 最后三年的图形:

plot(last(xts.gs, "3 years"), type="l", 
     main="US Federal Bonds 1 Year and 3 Years",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="auto",
     col=c("blue", "red"))
美国一年期与三年期国债利率周数据

图10.2: 美国一年期与三年期国债利率周数据

这些数据作为时间序列呈现出不平稳的表现。 对一年期利率作ACF估计:

acf(x1, main="")
一年期利率ACF

图10.3: 一年期利率ACF

这样的ACF是有缓慢变化趋势的典型表现。

对一年期利率序列作单位根检验:

ar(diff(x1), method="mle")
## 
## Call:
## ar(x = diff(x1), method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.3213   0.0372   0.0333   0.0577  -0.0093  -0.0008  -0.0813   0.0312  
##       9  
## -0.0469  
## 
## Order selected 9  sigma^2 estimated as  0.03103
fUnitRoots::adfTest(x1, lags=9, type="ct")
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 9
##   STATISTIC:
##     Dickey-Fuller: -2.3813
##   P VALUE:
##     0.4169 
## 
## Description:
##  Tue Mar 03 18:51:45 2020 by user: user

单位根检验结果不显著, 说明一年期利率序列是单位根过程。

如果以\(x_{3t}\)为因变量, 以\(x_{1t}\)为自变量作线性回归: \[ x_{3t} = \beta_0 + \beta_1 x_{1t} + e_t \] 就会遇到\(\{ e_t \}\)序列相关的问题, 甚至于\(\{ e_t \}\)可能有单位根造成方差线性增长, 有异方差问题。

先用散点图考察两个序列的关系:

plot(x1, x3, type="p", pch=16, cex=0.5,
     xlab="1 Year", ylab="3 Years",
     main="US Federal Bonds 3 Years Rate Vs. 1 Year Rate")
美国一年期与三年期国债利率周数据散点图

图10.4: 美国一年期与三年期国债利率周数据散点图

散点图显示同期的一年期利率和三年期利率高度相关。

试拟合线性回归模型:

lm1 <- lm(x3 ~ x1); summary(lm1)
## 
## Call:
## lm(formula = x3 ~ x1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.82319 -0.37691 -0.01462  0.38661  1.35679 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.83214    0.02417   34.43   <2e-16 ***
## x1           0.92955    0.00357  260.40   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5228 on 2465 degrees of freedom
## Multiple R-squared:  0.9649, Adjusted R-squared:  0.9649 
## F-statistic: 6.781e+04 on 1 and 2465 DF,  p-value: < 2.2e-16

注意:不懂统计的人会盲目相信统计软件给出的如此专业的输出结果。 \(R^2 = 0.9649\)是一个十分令人欣喜的回归结果。 但是, 只有学习了相应的统计模型的有关知识并且知道模型的局限, 才能正确地运用统计模型。 上面的系数估计、\(\sigma_e^2\)估计、SE估计、\(t\)检验等等都是基于随机误差项独立同分布假定的。 如果假定不成立,\(\sigma_e^2\)的估计会严重低估,SE估计不相合,检验不准确。 在经济和金融数据分析中一定要对回归模型结果进行诊断分析, 最常用的是残差分析, 其中一项是残差序列相关性的检验。

传统的残差序列相关性检验有Durbin-Watson检验, 在时间序列分析软件的支持下, 我们可以做ACF图, 进行Ljung-Box白噪声检验。

acf(lm1$residuals, main="")
三年期对一年期利率回归残差的ACF

图10.5: 三年期对一年期利率回归残差的ACF

残差的ACF表明结果仍有很强的自相关,甚至于可能有单位根。 白噪声检验:

Box.test(lm1$residuals, lag=9)
## 
##  Box-Pierce test
## 
## data:  lm1$residuals
## X-squared = 19303, df = 9, p-value < 2.2e-16

结果显著,拒绝残差序列为白噪声的零假设。

对残差序列作单位根检验:

ar(diff(lm1$residuals), method="mle")
## 
## Call:
## ar(x = diff(lm1$residuals), method = "mle")
## 
## Coefficients:
##       1        2        3        4        5        6        7        8  
##  0.2072  -0.0030  -0.0100   0.0911  -0.0641  -0.0642  -0.0475   0.0193  
##       9       10       11       12  
##  0.0163  -0.0509  -0.0174   0.0706  
## 
## Order selected 12  sigma^2 estimated as  0.005036
fUnitRoots::adfTest(lm1$residuals, lags=12, type="nc")
## Warning in fUnitRoots::adfTest(lm1$residuals, lags = 12, type = "nc"): p-value
## smaller than printed p-value
## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 12
##   STATISTIC:
##     Dickey-Fuller: -4.1201
##   P VALUE:
##     0.01 
## 
## Description:
##  Tue Mar 03 18:51:57 2020 by user: user

结果显著,拒绝残差有单位根的零假设。 尽管如此,鉴于残差的高度自相关, 还是不能轻易使用线性回归。

如果回归模型\(y_t = \beta_0 + \beta_1 x_t + e_t\)\(y_t\)\(x_t\)都是单位根过程而\(e_t\)是线性时间序列, 则称\(\{ y_t \}\)\(\{ x_t \}\)序列存在协整关系(cointegrated)。 在目前的数据中, 因为利率序列有单位根, 残差虽然没有单位根但是ACF衰减很慢, 不符合协整的要求。

为此, 对利率序列作一阶差分,考虑利率差分之间的关系。 分别记差分后的两个序列为\(\{c_{1t}\}\)\(\{c_{3t}\}\)

c1 <- diff(x1)
c3 <- diff(x3)

一年期利率的差分序列的ACF:

acf(c1, main="")
一年期利率差分序列的ACF

图10.6: 一年期利率差分序列的ACF

这个ACF的表现已经不是单位根过程, 符合线性时间序列的ACF的表现。

一年期与三年期的两个差分序列的散点图:

plot(c1, c3, type="p", pch=16, cex=0.5,
     xlab="1 Year", ylab="3 Years",
     main="US Federal Bonds 3 Years Rate Diff1 Vs. 1 Year Rate Diff1")
一年期利率差分序列的ACF

图10.7: 一年期利率差分序列的ACF

两者仍有强烈的线性相关。作线性回归:

lm2 <- lm(c3 ~ c1); summary(lm2)
## 
## Call:
## lm(formula = c3 ~ c1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42459 -0.03578 -0.00117  0.03467  0.48921 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0001051  0.0013890  -0.076     0.94    
## c1           0.7919323  0.0073391 107.906   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06897 on 2464 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8253 
## F-statistic: 1.164e+04 on 1 and 2464 DF,  p-value: < 2.2e-16

结果表面上也很好, 但仍需作序列相关性诊断分析。 残差的ACF:

acf(lm2$residuals, main="")
差分之间回归残差的ACF

图10.8: 差分之间回归残差的ACF

这个ACF表现得与线性时间序列ACF相像, 但是滞后1处仍显著不等于零。 作Ljung-Box白噪声检验:

Box.test(lm2$residuals, lag=9)
## 
##  Box-Pierce test
## 
## data:  lm2$residuals
## X-squared = 129.57, df = 9, p-value < 2.2e-16

结果显著, 说明回归模型中的误差项存在序列自相关, 这样, 上面的两个差分序列的回归结果中仅两个回归系数估计是可信的, 误差方差估计、SE估计、假设检验结果均不可用。

因为上面的回归残差的ACF主要在滞后1显著,其它位置虽然也有显著但值不大, 所以考虑\(e_t\)用MA(1), 回归变成

\[ c_{3t} = \beta_1 c_{1t} + e_t, \quad e_t = \varepsilon_t + \theta \varepsilon_{t-1} \] 其中\(\{ \varepsilon_t \}\)为零均值独立同分布白噪声列。 这里没有常数项是因为两个序列都是差分过的, 如果原来的关系有常数项会在差分过程中被消去。

R的arima()函数提供了将外生回归自变量与AMRA模型一同估计的功能。 用xreg=指定外生回归自变量。

mres <- arima(c3, xreg=c1, order=c(0,0,1), include.mean=FALSE)
mres
## 
## Call:
## arima(x = c3, order = c(0, 0, 1), xreg = c1, include.mean = FALSE)
## 
## Coefficients:
##          ma1      c1
##       0.1823  0.7936
## s.e.  0.0196  0.0075
## 
## sigma^2 estimated as 0.0046:  log likelihood = 3136.62,  aic = -6267.23

这拟合了如下的带有MA(1)误差的无截距项的一元线性回归模型:

\[ c_{3t} = 0.7936 c_{1t} + e_t, \quad e_t = \varepsilon_t + 0.1823 \varepsilon_{t-1}, \ \hat\sigma_{\varepsilon}^2 = 0.0046 \]

\(\hat\beta_1\)\(\hat\theta\)的SE估计来看, 两个系数都是显著不为零的。

小结:带时间序列误差的线性回归模型建模步骤

  1. 拟合一个线性回归模型,并检验残差的序列相关性
  2. 如果残差序列是单位根非平稳的,则对因变量和自变量都做一阶差分。 然后对差分后的序列再进行第一步。 如果残差序列是平稳的, 则对残差序列识别一个ARMA模型, 并相应地修改线性回归模型。
  3. 用最大似然估计法对回归模型与ARMA模型进行联合估计, 并对模型进行检验, 看是否需要改进。 主要可以使用Ljung-Box对残差进行白噪声检验。