北京大学金融时间序列分析讲义第15章： 线性时间序列案例学习「美国月失业率」

da <- read_table("m-unrate.txt", col_types = cols(.default = col_double()))
xts.unrate <- xts(
  da[["rate"]], with(da, make_date(Year, mon, dd)))
tclass(xts.unrate) <- "yearmon"
ts.unrate <- ts(da[["rate"]], start=c(1948,1), frequency = 12)
ts.dur <- diff(ts.unrate)

plot(
  xts.unrate, ylim=c(0, 12),
  main="US Monthly Unemployment Rate",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years",
  col="red"
)

图15.1: 美国月度失业率（经季节调整后）

plot(
  xts.unrate["1948/1968"], ylim=c(0, 12),
  main="US Monthly Unemployment Rate",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years",
  col="red"
)

图15.2: 美国月度失业率（经季节调整后）1948-1968

plot(
  xts.unrate["1968/1988"], ylim=c(0, 12),
  main="US Monthly Unemployment Rate",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years",
  col="red"
)

图15.3: 美国月度失业率（经季节调整后）1968-1988

plot(
  xts.unrate["1988/"], ylim=c(0, 12),
  main="US Monthly Unemployment Rate",
  major.ticks="years", minor.ticks=NULL, 
  grid.ticks.on="years",
  col="red"
)

图15.4: 美国月度失业率（经季节调整后）1988-2010

acf(ts.unrate, lag.max=36, main="")

图15.5: 美国月度失业率的ACF

forecast::Acf(ts.dur, lag.max=36, main="")

图15.6: 美国月度失业率差分的ACF

pacf(ts.dur, lag.max=36, main="")

图15.7: 美国月度失业率差分的PACF

fUnitRoots::adfTest(ts.unrate, lags=12, type="ct")

resm <- arima(
  ts.unrate, order=c(1,1,5),
  seasonal=list(order=c(1,0,1), period=12)
)
resm

mod1 <- arima(
  ts.unrate, order=c(1,1,5),
  seasonal=list(order=c(1,0,1), period=12),
  fixed=c(NA, NA, NA, 0, 0, NA, NA, NA)
)
mod1

图15.8: 失业率ARIMA模型的残差检验

resm2 <- arima(
  ts.unrate, order=c(0,1,0),
  seasonal=list(order=c(1,0,1), period=12)
)
resm2

forecast::Acf(resm2$residuals, lag.max=36, main="")

图15.9: 单变量仅季节项的残差ACF

pacf(resm2$residuals, lag.max=36, main="")

图15.10: 单变量仅季节项的残差PACF

resm3 <- arima(
  ts.unrate, order=c(5,1,0),
  seasonal=list(order=c(1,0,1), period=12)
)
resm3

fixed <- rep(NA_real_, 7); fixed[1] <- 0
mod2 <- arima(
  ts.unrate, order=c(5,1,0),
  seasonal=list(order=c(1,0,1), period=12),
  fixed=fixed
)

图15.11: 单变量候选模型检验

pred1 <- function(y=ts.unrate, n.pred=73){
  T <- length(y)
  ypred <- rep(NA_real_, T)
  for(h in (T-n.pred):(T-1)){
    mod <- arima(
      y[1:h], order=c(1,1,5),
      seasonal=list(order=c(1,0,1), period=12),
      fixed=c(NA, NA, NA, 0, 0, NA, NA, NA)
    )
    ypred[h+1] <- predict(mod, n.ahead=1, se.fit=FALSE)
  }
  
  ypred
}

pred2 <- function(y=ts.unrate, n.pred=73){
  T <- length(y)
  ypred <- rep(NA_real_, T)
  fixed <- rep(NA_real_, 7); fixed[1] <- 0
  for(h in (T-n.pred):(T-1)){
    mod <- arima(
      y[1:h], order=c(5,1,0),
      seasonal=list(order=c(1,0,1), period=12),
      fixed=fixed, transform.pars=FALSE
    )
    ypred[h+1] <- predict(mod, n.ahead=1, se.fit=FALSE)
  }
  
  ypred
}

pred.err <- function(y, n.pred=length(ypred), ypred){
  T <- length(y)
  i1 <- T - n.pred + 1
  RMFSE <- sqrt( mean((y[i1:T] - ypred[i1:T])^2) )
  MAFE <- mean(abs(y[i1:T] - ypred[i1:T]))
  
  c(RMFSE=RMFSE, MAFE=MAFE)
}

err1 <- pred.err(y=ts.unrate, n.pred=73, 
                 ypred=pred1(ts.unrate, n.pred=73))
err2 <- pred.err(y=ts.unrate, n.pred=73, 
                 ypred=pred2(ts.unrate, n.pred=73))
err.tab1 <- data.frame(
  Model=c("mod1", "mod2"),
  RMFSE=c(err1[["RMFSE"]], err2[["RMFSE"]]),
  MAFE=c(err1[["MAFE"]], err2[["MAFE"]])
)
knitr::kable(err.tab1, digits=4)

da2 <- read_table("m-unrateic.txt", 
                   col_types = cols(.default = col_double()))
da2[,5:9] <- da2[,5:9]/1000
xts.unic <- xts(
  da2[,-(1:3)], with(da2, make_date(year, mon, dd)))
tclass(xts.unic) <- "yearmon"
unrate <- da2[["rate"]]
head(da2)

plot(xts.unic[,c("rate", "icm1")],
     multi.panel=TRUE, yaxis.same=FALSE,
     main="月失业率与上个月申请失业救济金人数（单位：百万人）",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="years",
     col=c("red", "green"))

图15.12: 月失业率与上个月申请失业救济金人数（单位：百万人）

plot(da2[["icm1"]], da2[["rate"]], type="p",
     pch=16, col="gray",
     xlab="上个月申请失业救济金人数（单位：百万人）",
     ylab="月失业率")

图15.13: 月失业率对上个月申请失业救济金人数

cor(da2[["icm1"]], da2[["rate"]])

使用前一月申请人数总和

图15.13给出了从1967年2月到2010年9月的月失业率xt的时序图，
以及对应的前一个月的首次申请失业救济金人数的ct−1序列的时序图。
可以看出两者呈现一种共同的运动。

先拟合一个简单的一元线性回归模型：

lm1 <- lm(rate ~ icm1, data=da2)
summary(lm1)

## 
## Call:
## lm(formula = rate ~ icm1, data = da2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8638 -0.7008 -0.1299  0.7366  3.1746 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.5202     0.1785   8.518   <2e-16 ***
## icm1          2.9047     0.1097  26.475   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.051 on 522 degrees of freedom
## Multiple R-squared:  0.5731, Adjusted R-squared:  0.5723 
## F-statistic: 700.9 on 1 and 522 DF,  p-value: < 2.2e-16

结果的回归复相关系数平方为0.57，
说明有一定预报能力。
模型为

xt=1.52+2.905ct−1+nt

其中nt为误差项。

但是，这个回归有两个问题：

• 残差很可能有序列自相关，
  使得回归的检验是不可用的，
  估计也不是最有效的；
• 如果回归因变量和自变量存在单位根非平稳，
  结果可能是虚假的。

回归残差的ACF:

lm1 <- lm(rate ~ icm1, data=da2)
acf(lm1$residuals, lag.max=36, main="")

残差的ACF衰减缓慢，
可能由于原始序列有单位根。

回归残差的PACF:

lm1 <- lm(rate ~ icm1, data=da2)
pacf(lm1$residuals, lag.max=36, main="")

前3个PACF很较高，后面也有多处超出界限。
从回归残差的ACF和PACF来看不一定能用ARMA模型建模。

为克服模型设定上的困难，
我们应用试错法和一些先验知识。

1. 由于是月度数据，
   模型可能需要季节项。
2. 失业率呈现出一定的非固定的周期波动，
   来源于商业周期，
   这表明失业率模型的一些特征根是复根，
   因此AR的阶在2以上；
3. 回归残差的ACF都是正的，
   可考虑MA阶为整数，试选取q=3。

暂定失业率带有月申请人数解释变量的模型为：

=(1−ϕ1B−ϕ2B2)(1−ΦB12)(xt−β0−β1ct−1)(1+θ1B+θ2B2+θ3B3)(1+ΘB12)at

resm4 <- arima(
  da2[["rate"]], xreg=da2[["icm1"]], 
  order=c(2,0,3),
  seasonal=list(order=c(1,0,1), period=12)
); resm4

## 
## Call:
## arima(x = da2[["rate"]], order = c(2, 0, 3), seasonal = list(order = c(1, 0, 
##     1), period = 12), xreg = da2[["icm1"]])
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     ma3    sar1     sma1  intercept
##       1.8997  -0.9021  -0.8932  0.1458  0.0555  0.6501  -0.8520     6.0373
## s.e.  0.0332   0.0331   0.0543  0.0565  0.0466  0.0824   0.0586     0.3706
##       da2[["icm1"]]
##              0.0772
## s.e.         0.0212
## 
## sigma^2 estimated as 0.02419:  log likelihood = 227.7,  aic = -435.39

因为数据不同，这里的AIC不与前面的单变量模型比较。
结果表明ma3系数不显著，所以MA改为q=2建模：

mod3 <- arima(
  da2[["rate"]], xreg=da2[["icm1"]], 
  order=c(2,0,2),
  seasonal=list(order=c(1,0,1), period=12)
); mod3

## 
## Call:
## arima(x = da2[["rate"]], order = c(2, 0, 2), seasonal = list(order = c(1, 0, 
##     1), period = 12), xreg = da2[["icm1"]])
## 
## Coefficients:
##          ar1      ar2      ma1     ma2    sar1     sma1  intercept
##       1.9123  -0.9145  -0.9100  0.1860  0.6465  -0.8483     6.1111
## s.e.  0.0283   0.0282   0.0527  0.0479  0.0823   0.0591     0.3748
##       da2[["icm1"]]
##              0.0782
## s.e.         0.0213
## 
## sigma^2 estimated as 0.02426:  log likelihood = 226.97,  aic = -435.93

各个系数都显著。
模型的残差诊断：

除了残差有较多异常值以外，
残差呈现出白噪声性质。

模型mod3可以写成：

=(1−1.91B+0.91B2)(1−0.65B12)(xt−6.11−0.0782ct−1)(1−0.91B+0.18B2)(1−0.84B12)at, Var(at)=0.02426, AIC=−435.93(15.3)

模型中季节项的两个系数比较接近，
这往往可以产生部分长记忆的现象。

模型中AR部分除了有一个对应于周期12的特征根以外，
其它的特征根和对应的周期为：

polyroot(c(1, -coef(mod3)[1:2]))

## [1] 1.045522+0.019388i 1.045522-0.019388i

2*pi/Arg(polyroot(c(1, -coef(mod3)[1:2]))[1])/12

## [1] 28.23933

模型代表的周期为28个月。

使用前一月各周申请人数

这里使用前一个月的前四周的申请人数都作为自变量。
下图为失业率与前一个月的前两周的申请人数的数据，
三个序列走势是相近的：

plot(xts.unic[,c("rate", "w1m1", "w2m1")],
     multi.panel=TRUE, yaxis.same=FALSE,
     main="月失业率与上个月前两周申请失业救济金人数（单位：百万人）",
     major.ticks="years", minor.ticks=NULL, 
     grid.ticks.on="years",
     col=c("red", "green", "cyan"))

用四个解释变量作多元线性回归：

lm2 <- lm(rate ~ w1m1 + w2m1 + w3m1 + w4m1, data=da2)
summary(lm2)

## 
## Call:
## lm(formula = rate ~ w1m1 + w2m1 + w3m1 + w4m1, data = da2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.40629 -0.67602 -0.06686  0.62419  2.55888 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5160     0.1651   3.125 0.001877 ** 
## w1m1          6.5221     2.1145   3.084 0.002148 ** 
## w2m1          9.6711     2.8630   3.378 0.000785 ***
## w3m1         -2.4455     2.7980  -0.874 0.382506    
## w4m1          1.6626     2.0624   0.806 0.420528    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8841 on 519 degrees of freedom
## Multiple R-squared:  0.6998, Adjusted R-squared:  0.6975 
## F-statistic: 302.5 on 4 and 519 DF,  p-value: < 2.2e-16

复相关系数平方为0.70，
高于仅使用总数作为解释变量时的0.57。
当然，
线性回归每增加一个解释变量其复相关系数平方都会增加，
有过度拟合危险。
这个回归的可能问题与用总数作为解释变量的情形类似，
另外，
模型中后第三、四周的解释变量不显著。
这并不表明第三、四周的申请人数对下个月失业率没有预报能力，
而是其预报能力可以被第一、二周的数据代替。
比如，
仅使用第三、四周建模：

lm2b <- lm(rate ~ w3m1 + w4m1, data=da2)
summary(lm2b)

## 
## Call:
## lm(formula = rate ~ w3m1 + w4m1, data = da2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.62420 -0.69034 -0.07828  0.63761  2.89772 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.6694     0.1707   3.920   0.0001 ***
## w3m1         10.2363     2.1342   4.796 2.11e-06 ***
## w4m1          4.7578     2.1012   2.264   0.0240 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9226 on 521 degrees of freedom
## Multiple R-squared:  0.6718, Adjusted R-squared:  0.6706 
## F-statistic: 533.3 on 2 and 521 DF,  p-value: < 2.2e-16

结果仍是显著的。

下面仅使用第一、二周的申请人数预报下个月的失业率：

lm3 <- lm(rate ~ w1m1 + w2m1, data=da2)
summary(lm3)

## 
## Call:
## lm(formula = rate ~ w1m1 + w2m1, data = da2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.40653 -0.66933 -0.05949  0.63075  2.56954 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.5127     0.1647   3.113  0.00195 ** 
## w1m1          6.4594     2.0709   3.119  0.00191 ** 
## w2m1          8.9609     2.0872   4.293  2.1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8832 on 521 degrees of freedom
## Multiple R-squared:  0.6993, Adjusted R-squared:  0.6981 
## F-statistic: 605.8 on 2 and 521 DF,  p-value: < 2.2e-16

复相关系数平方仍为0.70，
相比使用四个解释变量的模型基本没有下降。

仅使用前两周的数据预报还有一个额外的好处，
就是给出预报的时间提前了，
可以在前一个月15号的时候就预报下个月的失业率，
如果使用总数则需要在月底才能预报下个月1号报告的失业率。

考察模型的残差ACF：

acf(lm3$residuals, lag.max=36, main="")

回归残差ACF衰减缓慢，
怀疑有单位根残余。
相关系数都为正数。

pacf(lm3$residuals, lag.max=36, main="")

这个PACF比较符合截尾特征，
取p=3可能合适。
使用与mod3类似的设定：

mod4 <- arima(
  da2[["rate"]], xreg=as.matrix(da2[,c("w1m1", "w2m1")]), 
  order=c(2,0,2),
  seasonal=list(order=c(1,0,1), period=12)
); mod4

## 
## Call:
## arima(x = da2[["rate"]], order = c(2, 0, 2), seasonal = list(order = c(1, 0, 
##     1), period = 12), xreg = as.matrix(da2[, c("w1m1", "w2m1")]))
## 
## Coefficients:
##          ar1      ar2      ma1     ma2    sar1     sma1  intercept    w1m1
##       1.9172  -0.9197  -0.9958  0.2532  0.6111  -0.7915     5.6555  0.4265
## s.e.  0.0269   0.0268   0.0563  0.0507  0.1119   0.0883     0.3912  0.2721
##         w2m1
##       0.9693
## s.e.  0.3206
## 
## sigma^2 estimated as 0.024:  log likelihood = 230.29,  aic = -440.59

模型AIC为−440.59，小于mod3的−435.93。
样本内比较mod4优于mod3。

mod4模型残差诊断：

残差诊断可以接受。

模型mod4公式为：

=(1−1.92B+0.92B2)(1−0.61B12)(xt−5.66−0.4265w1,t−1−0.9693w2,t−1)(1−1.00B+0.25B2)(1−0.79B12)at, Var(at)=0.024, AIC=−4440.59(15.4)

仔细比较mod3和mod4，
可以发现回归残差的ARMA模型的系数基本相同。

mod5 <- arima(
  da2[["rate"]],
  order=c(1,1,2),
  seasonal=list(order=c(1,0,1), period=12)
); mod5

图15.21: 时间区间缩短的一元时间序列的残差诊断

err.tab2 <- data.frame(
  "模型"=c("使用月申请人数", "使用两周申请人数", "不使用解释变量"),
  AIC=c(-435.93, -440.59, -426.17))
knitr::kable(err.tab2)

pred3 <- function(da=da2, n.pred=53){
  T <- nrow(da)
  ypred <- rep(NA_real_, T)
  for(h in (T-n.pred):(T-1)){
    mod <- arima(
      da[["rate"]][1:h], xreg=da[["icm1"]][1:h], 
      order=c(2,0,2),
      seasonal=list(order=c(1,0,1), period=12)
    )
    ypred[h+1] <- predict(mod, n.ahead=1, se.fit=FALSE,
                          newxreg=cbind(da[["icm1"]][h+1]))
  }
  
  ypred
}

pred4 <- function(da=da2, n.pred=53){
  T <- nrow(da)
  ypred <- rep(NA_real_, T)
  for(h in (T-n.pred):(T-1)){
    mod <- arima(
      da[["rate"]][1:h], xreg=as.matrix(da2[1:h,c("w1m1", "w2m1")]), 
      order=c(2,0,2),
      seasonal=list(order=c(1,0,1), period=12)
    )
    ypred[h+1] <- predict(mod, n.ahead=1, se.fit=FALSE,
                          newxreg=cbind(da[["w1m1"]][h+1], da[["w2m1"]][h+1]))
  }
  
  ypred
}

pred5 <- function(da=da2, n.pred=53){
  T <- nrow(da)
  ypred <- rep(NA_real_, T)
  for(h in (T-n.pred):(T-1)){
    mod <- arima(
      da[["rate"]][1:h],
      order=c(1,1,2),
      seasonal=list(order=c(1,0,1), period=12)
    )
    ypred[h+1] <- predict(mod, n.ahead=1, se.fit=FALSE)
  }
  
  ypred
}

pred.tab3 <- data.frame(
  "模型"=c("使用月申请人数", "使用两周申请人数", "不使用解释变量"),
  AIC=c(-435.93, -440.59, -426.17),
  RMFSE=rep(0.0, 3),
  MAFE=rep(0.0, 3)
)
funcs <- list(pred3, pred4, pred5)
for(i in 1:3){
  predf <- funcs[[i]]
  ypred <- predf(da=da2, n.pred=53)
  err <- pred.err(y=da2[["rate"]], n.pred=53, ypred=ypred)
  pred.tab3[i, "RMFSE"] <- err["RMFSE"]
  pred.tab3[i, "MAFE"] <- err["MAFE"]
}

knitr::kable(pred.tab3, digits=4)

abs(polyroot(c(1, -coef(mod3)[1:2])))

abs(polyroot(c(1, -coef(mod4)[1:2])))