layout: true --- class: inverse, center, middle background-image: url(figures/titlepage16-9.png) background-size: cover <br> # 商业预测 ## 第五讲:动态回归模型及应用 <img src="figures/qr.png" width="150px"/> #### 康雁飞 | 北航经管学院 | 2020年秋 --- # 核心思路 - 在 ARIMA 模型中,我们仅利用变量的实际观测值进行建模,却没有考虑到随时间的推移而产生的相关信息的影响。 - 在回归模型中,模型考虑了大量预测变量的相关信息,但并未采用ARIMA模型处理时间序列中的动态因素。 - 在本讲中,我们考虑如何扩展 ARIMA 模型,将一些与预测变量相关的信息纳入模型之中。 --- # 回归模型回顾 我们曾建立如下形式的回归模型: `$$y_t = \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \varepsilon_t,$$` 其中, `\(y_t\)` 是 `\(k\)` 个预测变量( `\(x_{1,t},\dots,x_{k,t}\)` )的线性组合, `\(\varepsilon_t\)` 是不相关的误差项(例如,误差项为白噪声)。 --- class: inverse, center, middle # 动态回归模型 --- # 动态回归模型 我们现在允许误差项中存在自相关。我们用 `\(\eta_t\)` 代替 `\(\varepsilon_t\)`。 `\(\varepsilon_t\)` 表示误差项,且利用 ARIMA 模型拟合残差 `\(\eta_t\)` 。例如,若残差 `\(\eta_t\)` 满足 ARIMA(1,1,1),则回归模型可写为: `\begin{align*} y_t &= \beta_0 + \beta_1 x_{1,t} + \dots + \beta_k x_{k,t} + \eta_t,\\ & (1-\phi_1B)(1-B)\eta_t = (1+\theta_1B) \varepsilon_t, \end{align*}` 其中 `\(\varepsilon_t\)` 为白噪声。模型中存在两个残差项:一是来自回归模型的误差,可用 `\(\eta_t\)` 表示;二是来自 ARIMA 模型的误差,可用 `\(\varepsilon_t\)` 表示。且假设来自 ARIMA 模型的误差为白噪声。 --- # 例:美国个人消费及收入 .pull-left[ ```r cbind("消费变化" = uschange[, 1], "收入变化" = uschange[, 2]) %>% autoplot(facets=TRUE) + xlab("年份") + ylab("") + ggtitle("美国个人消费和收入的季度变化")+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/usconsump-out-1.png" style="display: block; margin: auto;" /> ] --- # 动态回归模型的拟合 ```r (fit <- auto.arima(uschange[,"Consumption"], xreg=uschange[,"Income"])) ``` ``` ## Series: uschange[, "Consumption"] ## Regression with ARIMA(1,0,2) errors ## ## Coefficients: ## ar1 ma1 ma2 intercept xreg ## 0.6922 -0.5758 0.1984 0.5990 0.2028 ## s.e. 0.1159 0.1301 0.0756 0.0884 0.0461 ## ## sigma^2 estimated as 0.3219: log likelihood=-156.95 ## AIC=325.91 AICc=326.37 BIC=345.29 ``` --- # 动态回归模型的拟合 拟合的模型为: `$$\begin{array}{l}{y_{t}=0.60+0.20 x_{t}+\eta_{t}} \\ {\eta_{t}=0.69 \eta_{t-1}+\varepsilon_{t}-0.58 \varepsilon_{t-1}+0.20 \varepsilon_{t-2}} \\ {\varepsilon_{t} \sim \mathrm{NID}(0,0.322)}\end{array}$$` --- # 动态回归模型的误差 .pull-left[ ```r cbind("回归误差" = residuals(fit, type="regression"), "ARIMA误差" = residuals(fit, type="innovation")) %>% autoplot(facets=TRUE)+ xlab('年份')+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/usconsumpres-out-1.png" style="display: block; margin: auto;" /> ] --- # 动态回归模型的误差 <div class="figure" style="text-align: center"> <img src="BF-L5-dynamic_files/figure-html/usconsumpres2-1.png" alt="误差(即ARIMA模型误差)与白噪声没有明显差异。" /> <p class="caption">误差(即ARIMA模型误差)与白噪声没有明显差异。</p> </div> ``` ## ## Ljung-Box test ## ## data: Residuals from Regression with ARIMA(1,0,2) errors ## Q* = 5.8916, df = 3, p-value = 0.117 ## ## Model df: 5. Total lags used: 8 ``` --- # 预测 .pull-left[ 我们将会预测未来八个季度的预测值,未来的个人可支配收入的百分比变化等于过去四十年间个人可支配收入的平均百分比变化。 ```r fcast <- forecast(fit, xreg=rep(mean(uschange[,2]),8)) autoplot(fcast) + xlab("年份") + ylab("百分比变化")+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/usconsump3-out-1.png" style="display: block; margin: auto;" /> ] --- # 例:耗电量的预测 <div class="figure" style="text-align: center"> <img src="BF-L5-dynamic_files/figure-html/elecscatter-1.png" alt="2014年澳大利亚维多利亚州每日耗电量与每日最高温度的散点图。" /> <p class="caption">2014年澳大利亚维多利亚州每日耗电量与每日最高温度的散点图。</p> </div> --- # 耗电量与气温 .pull-left[ 数据集`elecdaily`中包含日度耗电量、工作日的指示变量(1表示为工作日,0表示非工作日)和当日最高温度。由于数据存在周度季节性,因此指定频率为7。 ```r cbind("耗电量" = elecdaily[,"Demand"], "温度" = elecdaily[,"Temperature"]) %>% autoplot(facets=TRUE)+ ylab("") + guides(colour="none")+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/electime-out-1.png" style="display: block; margin: auto;" /> ] --- # 建模 ```r xreg <- cbind(MaxTemp = elecdaily[, "Temperature"], MaxTempSq = elecdaily[, "Temperature"]^2, Workday = elecdaily[, "WorkDay"]) fit <- auto.arima(elecdaily[, "Demand"], xreg = xreg) ``` --- # 预测 .pull-left[ ```r fcast <- forecast(fit, xreg = cbind(rep(26,14), rep(26^2,14), c(0,1,0,0,1,1,1,1,1,0,0,1,1,1))) autoplot(fcast) + ylab("日度耗电量(千兆瓦)") + theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/elecdailyf-out-1.png" style="display: block; margin: auto;" /> ] --- class: inverse, center, middle # 滞后预测变量 --- # 滞后预测变量 有些时候,预测变量对待预测变量的影响不是简单而迅速的。例如,投放广告之后一段时间才会影响销售接受,而当月的销售收入会取决于过去几个月的广告支出。类似的,当工厂的安全政策的变化会立即减少工厂事故的发生,但随着时间的推移,工人会降低对新出安全政策的重视程度。 在这些情况下,我们需要再模型汇总引入预测变量的滞后项。假如模型中只有一个预测变量,则假如该预测变量的滞后项之后,滞后模型可写为: `$$y_t = \beta_0 + \gamma_0x_t + \gamma_1 x_{t-1} + \dots + \gamma_k x_{t-k} + \eta_t,$$` 其中, `\(\eta_t\)` 为 ARIMA 过程。并通过AICc中来选择 `\(K\)` 值和ARIMA过程中的 `\(p\)` 和 `\(q\)`。 --- # 例:电视广告和投保数量 .pull-left[ 美国一家保险公司希望通过在国家电视台投放广告来增加保险的销售量(且通常为新出的保险种类)。 ```r cbind("保险销量" = insurance[, "Quotes"], "广告支出" = insurance[, "TV.advert"]) %>% autoplot(facets=TRUE)+ xlab("年份") + ylab("") + ggtitle("保险销量和广告支出")+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/tvadvert-out-1.png" style="display: block; margin: auto;" /> ] --- # 滞后模型 在滞后模型中,我们考虑当月的广告和之前三个月的广告对报销销量的影响。并采用不存在滞后项的模型与滞后模型进行对比。 ```r # Lagged predictors. Test 0, 1, 2 or 3 lags. Advert <- cbind( AdLag0 = insurance[,"TV.advert"], AdLag1 = stats::lag(insurance[,"TV.advert"],-1), AdLag2 = stats::lag(insurance[,"TV.advert"],-2), AdLag3 = stats::lag(insurance[,"TV.advert"],-3))[1:NROW(insurance),] # Restrict data so models use same fitting period fit1 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1], d=0) fit2 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:2], d=0) fit3 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:3], d=0) fit4 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:4], d=0) ``` --- # 滞后模型 滞后我们通过各个模型的AICc值来选择预测变量的最优滞后阶数。 ```r c(fit1$aicc,fit2$aicc,fit3$aicc,fit4$aicc) ``` ``` ## [1] 68.49968 60.02357 62.83253 65.45747 ``` 由上述结果可以看出,当预测变量的滞后阶数为一时,AICc 值最小,因此最优模型中应当包含预测变量的一阶滞后。也就是说,应该采用当月和上月的广告支出对当月的保险销量进行建模。模型如下所示: ```r (fit <- auto.arima(insurance[,1], xreg=Advert[,1:2], d=0)) ``` ``` ## Series: insurance[, 1] ## Regression with ARIMA(3,0,0) errors ## ## Coefficients: ## ar1 ar2 ar3 intercept AdLag0 AdLag1 ## 1.4117 -0.9317 0.3591 2.0393 1.2564 0.1625 ## s.e. 0.1698 0.2545 0.1592 0.9931 0.0667 0.0591 ## ## sigma^2 estimated as 0.2165: log likelihood=-23.89 ## AIC=61.78 AICc=65.4 BIC=73.43 ``` --- # 滞后模型 模型的误差项为 AR(3)过程。模型可写为: `$$y_{t}=2.04+1.26 x_{t}+0.16 x_{t-1}+\eta_{t}$$` 其中, `\(y_t\)` 为第 `\(t\)` 月的保险销量, `\(x_t\)` 为第 `\(t\)` 月的广告支出, `$$\eta_{t}=1.41 \eta_{t-1}-0.93 \eta_{t-2}+0.36 \eta_{t-3}+\varepsilon_{t}$$` 其中, `\(\varepsilon_t\)` 为白噪声。 --- # 保险销量预测 .pull-left[假若我们已经知道未来某月的广告支出,那么我们可以预测出该月得保险销量。假设未来每月广告支出都为8,那么我们可以得到以下预测: ```r fc8 <- forecast(fit, h=20, xreg=cbind(AdLag0=rep(8,20), AdLag1=c(Advert[40,1],rep(8,19)))) autoplot(fc8) + ylab("保险销量") + ggtitle("当广告支出为8个单位时保险销量的预测值")+ theme(plot.title = element_text(hjust = 0.5)) ``` ] .pull-right[ <img src="BF-L5-dynamic_files/figure-html/tvadvertf8-out-1.png" style="display: block; margin: auto;" /> ] --- # 实验5 数据集`advert`中包含汽车零部件的月度销售额和广告支出。 1. 使用`autoplot`函数绘制上述数据。参数`facets=TRUE`可以起到什么样的作用? 2. 利用`tslm()`函数拟合标准回归模型 `\(y_t = a + b x_t + \eta_t\)`,其中 `\(y_t\)` 代表销售额, `\(x_t\)` 代表广告支出。 3. 证明残差具有显著的自相关性。 4. 使用`auto.arima()`函数重新拟合。误差模型的参数估计有什么差异?应选用哪种 ARIMA 模型拟合残差? 5. 检验拟合模型的残差。 6. 假设未来六个月的广告支出为每月十个单位,通过模型预测月度销售额,并生成预测区间。 --- # 总结 1. 动态回归模型 2. 滞后模型 3. 适用场景